function P = ipdivdiff(x,xdata,fdata)

%
% Interpolate Using Newton's Divided Difference (divdiff.m)
%
F = diag(divdiff(xdata,fdata))
n = length(xdata);
P = F(1)*ones(size(x));

for i = 2:n
  p = 1;
  for j = 1:(i-1)
    p = p .* (x-xdata(j));
  end
  p = F(i)*p;
  P = P + p
end