%Computes coefficients of interpolating polynomial using
%Newton's Divided differences.
%Input: x and y are vectors containing the x and y coordinates
%Output: coefficients c of interpolating polynomial in nested form
%Use with nest.m to evaluate interpolating polynomial
%
% Example of useage (with nest.m):
%  xp=[-1 0 2 3];  %x-values for interpolation
%  yp=[-5 -1 1 11]; %y-values for interpolation
%  [v,c]=newtdd(xp,yp); %c is the vector of coeffs, v is whole table
%  xx=linspace(-1,3);   %New points to evaluate Newton polynomial
%  yy=nest(c,xx,xp);
%  plot(xp,yp,'r*',xx,yy);  %Plot data points and Newton poly


function [v,c]=newtdd(x,y)
n=length(x);
for j=1:n
  v(j,1)=y(j);        % Fill in y column of Newton triangle
end
for i=2:n             % For column i,
  for j=1:n+1-i       % fill in column from top to bottom
    v(j,i)=(v(j+1,i-1)-v(j,i-1))/(x(j+i-1)-x(j));
  end
end
for i=1:n
  c(i)=v(1,i);        % Read along top of triangle
end                   % for output coefficients