function  y = divdif(x,y)
%DIVDIFF Compute divided difference vector
%
%       d = DIVDIF(x,y) returns a vector of divided differences
%       associated with the points x and values y.  x and y must be
%       vectors of the same length and x must have distinct entries.
%       The returned vector d is the sequence of divided differences of
%       order 0, 1, ...
%
%           y(1), (y(2)-y(1))/(x(2)-x(1)), ...
%
%       These divided differences are the coefficient of Newton's
%       form of the Lagrange interpolating polynomial, interpolating
%       the values y(i) and the points x(i).  The function NINTERP to
%       evaluate this polynomial calls DIVDIFF.
%
%       See also NINTERP.

n = length(x);
if n ~= length(y) 
  error('divdif:  sizes of input vectors must agree')
end
for i = 1:n-1
  for j=n:-1:i+1
    y(j) = (y(j) - y(j-1)) / (x(j) - x(j-i));
  end
end

