function [P,A]=inth_pol(FUN,dFUN,xn,x);
%%  function [P,A]=inth_pol(FUN,dFUN,xn,x);
%%
%%  THE POLYNOMIAL INTERPOLATING the values and derivative values at xn
%%      Procedure: generate the divided difference table, determine the 
%%                 coefficients of the polynomial in Newton form and evaluate
%%                 using a nested Horner type algorithm
%%
%%  INPUT:
%%    xn   - the interpolation points as a row vector e.g. [0,1,2,3,4]
%%    FUN  - the function value 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
%%   dFUN  - derivative value data for the interpolation which may be entered 
%%           EITHER as a string, enclosed in ' ', with the name the m-file of 
%%           the derivative 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=ddd_tabl(FUN,dFUN,xn);  % get the divided difference table
n=length(xn)-1; % find the number of points
yn(1:2:2*n+1)=xn;
yn(2:2:2*n+2)=xn; % double the points
n=2*n+1;

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-yn(j))+A(j)*ones(size(x));
end;
