function yy = neville(xx,n,x,y)

 
% inputs:
%    n = order of interpolation (n+1 = # of points)
%    x(1),...,x(n+1)    x coords
%    y(1),...,y(n+1)    y coords
%    xx=evaluation point for interpolating polynomial p
%
% output:  yy=p(xx)

for i = n:-1:1
   for j = 1:i
      y(j) = (xx-x(j))*y(j+1) - (xx-x(j+n+1-i))*y(j);
      y(j) = y(j)/(x(j+n+1-i)-x(j));
   end
end

yy = y(1);