function [ v ] = newtoninterp( x,y,u )
%NEWTONINTERP Polynomial interpolation
%   Uses divided differences technique to get the Newton form
    v = y(1)*ones(size(u));
    p = ones(size(u));  % p = (u-x0)(u-x1)...(u-xi)
    for k=1:length(x)-1
        dx = x(k+1:end) - x(1:end-k);  % n-k vector of x(i+k)-x(i)
        y = diff(y)./dx;   % next row of divided difference table
        p = p.*(u-x(k));
        v = v + y(1)*p;
    end
end