function d = divdiff(x,y)
% d = divdiff(x,y)
% compute coefficients of Newton form of interpolating polynomial
% x: vector of nodes
% y: vector of y-values at nodes
% d: vector of Newton coefficients
%    d = [ f[x_1,...,x_n], f[x2,...,x_n], ..., f[x_n] ]
% use evnewt to evaluate interpolating polynomial
%
% Example: 
%   x = [0,1,2,4]; y = [1,2,3,1]; % given x and y values
%   d = divdiff(x,y)      % find coefficients of Newton form
%   xp = -1:.01:5;        % x-values for plotting
%   yp = evnewt(d,x,xp);  % evaluate Newton form at points xp
%   plot(x,y,'o',xp,yp)   % plot given pts and interpolating polynomial

n = length(x);
d = y;
for k=1:n-1
  for i=1:n-k
    d(i) = (d(i+1)-d(i))/(x(i+k)-x(i));
  end
end
