function d = divdif(x,f)
%
%  function d = divdif(x,f)
%
%  returns the vector of divided differences for the
%  Newton form of the interpolating polynomial
%
%  on entry:
%    x       abscisses, given as a column vector;
%    f       ordinates, given as a column vector.
%
%  on return:
%    d       divided differences
%
%  note:
%    MATLAB does not support negative or zero indices,
%    so this algorithm differs from the one on paper.
%
%  example (page 159 in Gerald-Wheatley):
%    x = [3.2 2.7 1.0 4.8 5.6]';
%    f = [22.0 17.8 14.2 38.3 51.7]';
%    d = divdif(x,f)
%
%  example (page 161 in Gerald-Wheatley):
%    x = [0.3 1.0 0.7 0.6 1.9 2.1]';
%    f = [-0.736 1.0 -0.104 -0.328 11.008 15.2120]';
%    d = divdif(x,f)
%
n = size(x,1);
d = f;
for i=2:n
  for j=1:i-1
    d(i) = (d(j) - d(i))/(x(j) - x(i));
  end;
end;