%USAGE D=divdiff(x,y)
% Input: 
%         x=(x1,...,xn)'
%         y=(y1,...,yn)'
% Output: 
%         D 
%For AMATH 3132 by Dr. J. Alam

function [D]=divdiff(x,y)
  % size of the tabulated data set
  n=length(x);
  % allocate memory for output matrix
  D=zeros(n,n+1);
  % calculate divided difference
  D(:,1)=x; D(:,2)=y;
  for k=3:n+1;
    for j=1:n+2-k;
      D(j,k)=(D(j+1,k-1)-D(j,k-1))/(D(j+k-2,1)-D(j,1));
    end
  end