%
% Compute Newton's interpolating polynomial
%
% function [a] = NewtIntPol( x, f, iprint )
%
% input:
%        x:     vector containing the 
%               interpolation points
%        f:     vector containing the values 
%               to be interpolated
%   iprint:     print flag
%               iprint = 0   no printout
%               iprint > 0   print divided difference table
%
% output:
%        a:     vector of coefficients of Newton's 
%               interpolating polynomial
%
%

function [a] = NewtIntPol( x, f, iprint );


n = size(x(:),1);


if size(f(:),1) ~= n 
   ifail = 1;
   return
end

A = zeros(n,n);

A(1:n,1) = f(:);
a(1)   = A(1,1);

for j = 2:n
    for i = j:n
        A(i,j) =  ( A(i,j-1) - A(i-1,j-1) ) / (x(i) - x(i-j+1)) ;
    end
    a(j) = A(j,j);
end

if( iprint > 0 )
    disp('        x         f  ')
    disp([x(:) A ])
end