%efficient algorithm for interpolation of points (x(i),y(i)), i=1,...,n
%fewer arithmetic operations compared to naive algorithm

n = 4;             %number of points
x = [-1;0;1;2];    %x,y are column vectors of length n
y = [-6;1;2;3];

a = zeros(n,n);    %this is array of divided differences
a(:,1) = y;        %first column of a is y(1),...,y(n)

for j=2:n
    for i=1:n-j+1
        a(i,j) = (a(i+1,j-1)-a(i,j-1))/(x(i+j-1)-x(i));
    end
end

disp('the coefficients are...')
for i=1:n
    a(1,i)         %divided differences f[x_0,...,x_i] for i=1,...,n
end