function [P,A]=int_poly(FUN,xn,x);
%%  function [P,A]=int_poly(FUN,xn,x);
%%
%%  THE POLYNOMIAL INTERPOLATING FUN at xn, generated from the divided
%%  difference table in Newton form and evaluated using a nested Horner 
%%  type algorithm
%%
%%  INPUT:
%%    xn   - the interpolation points as a row vector e.g. [0,1,2,3,4]
%%    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   - the points at which the interpolating polynomial is to be
%%           evaluated as a row vector e.g., x=0:.01:1
%%  OUTPUT
%%     P   - values of the polynomial at x
%%     A   - sequence of coefficients for the polynomial in Newton form
%%  

DD=dd_table(FUN,xn);  % get the divided difference table
n=length(xn)-1; % find the number of points

A=diag(DD(:,2:n+2)); % Get the coefficients of P from the top diagonal in
                    % the divided difference table

%%Now to evaluate the Newton polynomial with the Horner algorithm
P=A(n+1)*ones(size(x));
for j=n:-1:1;
      P=P.*(x-xn(j))+A(j)*ones(size(x));
end;
