function y = divdiff(x,y)
% in situ recursive computation of divided 
% differences for nodes x(i) and data y(i)
n = length(y)-1;
if (n > 0)
  y(1:n) = divdiff(x(1:n),y(1:n));
  for j=n:-1:1
    y(j) = (y(n+1)-y(j))/(x(n+1)-x(j));
  end
end

