function DD=DDtable(x,y)
n=length(x);
DD=zeros(n,n);
DD(:,1)=y(:);

for j=2:n
  for i=j:n
    DD(i,j) = (DD(i,j-1)-DD(i-1,j-1))/(x(i)-x(i-j+1));
  end
end