%algorithm for Newton interpolation of n points, avoiding extra storage
%first we compute coefficients, then we evaluate at user-specified point

n = 4;              %number of points

%first we compute coefficients of Newton interpolant

x = [-1,0,1,2];     %x = x-values for interpolation
y = [-6,1,2,3];     %y = y-values for interpolation

for j=2:n
    for i=n:-1:j
        y(i) = (y(i)-y(i-1))/(x(i)-x(i-j+1)); %old coefficients overwritten
    end
end                 %after this loop y = coefficients of Newton interpolant

disp('the coefficients are...')
y

%now we evaluate the interpolant at a specified point z

z = input('evaluate the interpolant at z = ')

p = 1;
w = y(1);

for i=2:n
    p = p*(z-x(i-1));
    w = w + y(i)*p;
end

disp('value of Newton interpolant at z is...')
w