function y = nev(xx,n,x,Q)
% NEV Neville's algorithm as a function: y= nev(xx,n,x,Q)
% where
%    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
%    y=interpolated/extrapolated value at xx
% (Ed Bueler; 2/4/02)

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);