% Newton's Divided Difference (forward) Interpolation
% Input: x=(x1,...,xn)'- x coordination vector
%        y=(y1,...,yn)'- y coordination vector
% Output: D is nx(n+1) upper triangular matrix elements are divided
% differences and it plots P(x) on interval [a,b].

function[D f]=divdiff(x,y)
n=length(x);
D=zeros(n,n+1);
D(:,1)=x; D(:,2)=y;
for k=3:n+1;
    for j=1:n+2-k;
        D(j,k)=(D(j+1,k-1)-D(j,k-1))/(D(j+k-2,1)-D(j,1));
    end
end
a=min(x)-2; b=max(x)+2; m=1000;
h=(b-a)/m;
for p=1:m;
    x=a+p*h; s=1;
    for k=n:-1:2
        s=s*D(1,k+1)*(x-D(k-1,1))+D(1,k);
    end
    f(p)=s;
end
plot(f)
___________________________________________________________________________
clear;
x=[0;2;3]
y=[1;3;6]
[D f]=divdiff(x,y);
 D

D =

         0    1.0000    1.0000    0.6667
    2.0000    3.0000    3.0000         0
    3.0000    6.0000         0         0