function dvided_diff
nodes=[-2,-1,0,1,2];  % nodes  
n=size(nodes);
n=n(2);          % determine the number of nodes 
diff=zeros(n,n); % set up a storage for the divided differences
diff(:,1)=(f(nodes))'; % sets up values of the function at the nodes = zeroth difference
for j=2:1:n
for i=1:1:n-j+1
diff(i,j)=(diff(i+1,j-1)-diff(i,j-1))/(nodes(i+j-1)-nodes(i));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the function versus interpolating polynomial 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(1,1,1);
clf; % clear the plot
% Plot the exact function 
a=min(nodes); %% figure the left end of the interval 
b=max(nodes); %% figure the right end of the interval 
x=[a:.1:b];
plot(x, f(x), '-r');
hold on; % to keep the plot on 
% calculate the approximated function 
m=size(x);
m=m(2);
y = zeros(1,m);
for i=1:1:m
y(i)=Poly(diff(1,:),nodes,x(i));
end
plot(x,y,'-b');
hold off;


function P = Poly(coeff, nodes, x)
n=size(nodes);
n=n(2);          % determine the number of nodes 
%it is expected that size of nodes is the same as the size of coeff! %
P=0;
for i=1:1:n
    q=1;
    for j=1:1:i-1
    q=q*(x-nodes(j));
    end
    P=P+coeff(i)*q;
end

function y=f(x);
y=sin(x);
%y=3.^x