function p_eval = interp ( x_nodes, divdif_y, x_eval )
%
% This is a function
%         p_eval = interp(x_nodes,divdif_y,x_eval)
% It calculates the Newton divided difference form of 
% the interpolation polynomial of degree m-1, where the 
% nodes are given in x_nodes, m is the length of x_nodes, 
% and the divided differences are given in divdif_y.  The 
% points at which the interpolation is to be carried out 
% are given in x_eval; and on exit, p_eval contains the 
% corresponding values of the interpolation polynomial.
%
%
%  Reference:
%
%    Kendall Atkinson, Weimin Han,
%    Elementary Numerical Analysis,
%    Wiley, 2004,
%    ISBN: 0471433373,
%    LC: QA297.A83.2004.
%
  m = length(x_nodes);
  p_eval = divdif_y(m)*ones(size(x_eval));
  for i=m-1:-1:1
    p_eval = divdif_y(i) + (x_eval - x_nodes(i)).*p_eval;
  end

  return
end

