function pn = polint ( Y, F )

%*****************************************************************************80
%
%MATLAB version of the fortran program on pages 68-69 of conte/deboor
% except that reading input and printing output is as well as the
% evaluation of the resulting polynomial in Newton form is left to the
% program that calls on POLINT.  

% Constructs of the Newton form  pn  for the polynomial of degree <= n
% that agrees with the data  F  about  f  at  Y = (Y(1), ..., Y(n+1)).
%
% Some or all of the interpolation sites may coincide, subject only
% to the following restrictions:
%   (1) If Y(i) == Y(i+k), then  Y(i) == Y(i+1) == ... == Y(i+k)
%   (2) If also Y(i-1) ~= Y(i) or if  i == 1, then
%       F(i+j) is the value of the jth derivative of  f  at Y(i),
%                                                          j=0:k . 
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
                 % initialize the Newton form structure 
   Y = Y(:); F = F(:); NP1 = length(Y); pn = struct('sites',Y,'coefs',F); 
   
                 % check the input sizes for compatibility
   if NP1 ~= length(F), disp('input arrays are inconsistent'), return, end
      
                 % construct divided differences
   for K=1:NP1-1
      FLAST = F(1);
      for I=1:NP1-K
         DY = Y(I+K)-Y(I);
         if (DY==0), F(I) = F(I+1)/K;
         else, F(I) = (F(I+1)-FLAST)/DY; FLAST = F(I+1);
         end
      end
      F(NP1-K+1) = FLAST;
   end
   pn.coefs = F;
