function interpolation_error
% This 
x_nodes = 0:0.00004:0.0002
y_values = fcn(x_nodes);
t = 0:.000001:max(x_nodes);

format long

for i = 2:length(x_nodes)
    
divdif_y = divdif(x_nodes(1:i),y_values(1:i));
p_eval = interp(x_nodes(1:i),divdif_y,t);


figure(1)
%axis([1 25 1 6])
hold all
%plot(t,ones(length(t)),'Linewidth',2,'Color','yellow')
plot(t,p_eval,'Linewidth',1,'Color','red')
%plot(x_nodes,y_values,'*','Color','black')
plot(t,fcn(t)*1.0001,'Linewidth',1,'Color','green')
end

%divdif_y(1,:)
fprintf('Newton divided differences are divdif=%g\n',divdif_y);
divdif_y(1,length(x_nodes))
fcn(x_nodes(length(x_nodes)))

figure(2)
title('propagation of rounding errors on equally distanced nodes')
plot(t,fcn(t)-p_eval,'Linewidth',2,'Color','green')

function fval = fcn(x)
fval = 1+ sin(3*x);
end
end