function P = neville(x,X,Y)
% function P = neville(x,X,Y)
% x = value approx function f(x) at (scalar)
% X = x-value nodes (row vector)
% Y = f(X) i.e. f sampled at the nodes (column vector)

n = length(X);
P = zeros(n,n);
P(:,1) = Y';

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

end 