function [d,D]=NewtonCoefficients(x,y)
% NEWTONCOEFFICIENTS divided differences for the Newton interpolation
%   [d,D]=NewtonCoefficients(x,y); computes the divided differences
%   needed for constructing the Newton form of the interpolating
%   polynomial through the points (x,y)

n=length(x)-1;                    % degree of interpolating polynomial
for i=1:n+1                       % divided differences
  D(i,1)=y(i);
  for j=1:i-1
    D(i,j+1)=(D(i,j)-D(i-1,j))/(x(i)-x(i-j));
  end
end
d=diag(D);
