function [Table,Px]=neville(FUN,xn,x);
%%function [Table,Px]=neville(FUN,xn,x);
%%
%% NEVILLE'S ALGORITHM for evaluating the interpolatory
%% polynomial P at x0,x1,...,xn for the values of FUN at the point x.
%%
%%  INPUT:
%%     xn    - n+1 points at which FUN is to be interpolated as a row
%%             vector  e.g. [1,2,3,4,5]
%%    FUN  - the data to be interpolated which may be entered EITHER as
%%           a string, enclosed in ' ', with the name the m-file of the 
%%           function ,  OR as a row vector of data in the same order as xn
%%     x     - point at which the interpolating polynomial is to be
%%             evaluated.
%%  OUTPUT:
%%     Table - matrix with first column xn, second column FUN(xn), and bottom 
%%             right hand entry equal to P(x).
%%     Px    - value of P(x) without the table.
%%  
%%

n=length(xn)-1;% find the number of interpolation points
Table(:,1)=xn';% make the first column equal to the given interpolation points

if isstr(FUN)==1;
   Table(:,2)=feval(FUN, xn');% 2nd column equal to function at xn
else
   Table(:,2)=FUN';
end;

for i=2:n+1;      %% The entries Q(i,j) will be Table(i+2,j+3)
   for j=3:i+1;   
         Table(i,j)=((x-xn(i))*Table(i-1,j-1)-(x-xn(i-j+2))*Table(i,j-1))...
                     /(xn(i-j+2)-xn(i));
   end;
end;
Px=Table(n+1,n+2);

