function y = nev(xx,n,x,Q)
% NEV   Neville's algorithm as a function
%       y= nev(xx,n,x,Q)
% inputs:
%    n = order of interpolation (n+1 = # of points)
%    x(1),...,x(n+1)    x coords
%    Q(1),...,Q(n+1)    y coords
%    xx=evaluation point for interpolating polynomial p
% output:  p(xx)
%
% example:
%    >> x = [0 0.1 0.4 0.5];      % n=3 case with n+1=4 points
%    >> y = cos(x);
%    >> P = nev(0.187654321,3,x,y) % value of polynomial P(xx)
%    >> f = cos(0.187654321)       % value of function f(xx)

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

y = Q(1);
