% Interpolation with Newton  polynomials

% Input: x,y: y=f(x)
%            xt: where interpolant is evaluated.
% Output: yt=f(xt)

function yt =  Newton_interpol(x,y,xt)

n = length(y); if length(x)~=n, error('x and y are not compatible'); end

c=y(:);
for j=2:n
    for i=n:-1:j
        c(i)=( c(i)-c(i-1) )  /  ( x(i) - x(i-j+1) );  % note that x(i)-x(i-1) is not right.
    end
end

% Nested evaluation of the polynomial
yt = c(n);
for i= n-1 :-1 :1
    yt = yt.*(xt - x(i)) + c(i);
end