function[d]=divdif(x,f)
%
% Computes divided difference table. x is the vector of interpolation
% points, f the vector of data.
%
n=length(x)-1;
d=f;
for i=2:n+1
  for j=n+1:-1:i
d(j)=(d(j)-d(j-1))/(x(j)-x(j-i+1));
  end
end


function p=interp(x,d,t)
%
% Computes interpolating polynomial at points in the vector t
% where x is the vector of nodes and d is the vector of divided
% differences produced by divdif.
%
n=length(x)-1;
k=length(t);
p=d(n+1)*ones(1,k);
for i=n:-1:1
p=d(i)+(t-x(i)).*p;
end
>> x=[0 1 2 3 4]

x =

     0     1     2     3     4

>> y=[2 4 8 10 6]

y =

     2     4     8    10     6

>> d= divdif(x,y)

d =

    2.0000    2.0000    1.0000   -0.6667         0

>> z=[1 1.5 2 2.5 3 3.5 4]

z =

    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000

>> p=interp(x,d,z)

p =

    4.0000    6.0000    8.0000    9.5000   10.0000    9.0000    6.0000