function y = intpolyvalnewt (x, xi)
%
% function y = intpolyvalnewt (x, xi)
%
% Valuta il polinomio interpolatore dei punti [xi(:,1),xi(:,2)] sui
% punti x, mediante la forma di Newton.
yi = xi(:, 2);
xi = xi(:, 1);
n = length (xi) - 1;
d = yi;
w = ones (size (x));
y = zeros (size (x));
for i = 1:n + 1
  for j = 1:i - 1
    d(i) = (d(i) - d(j)) / (xi(i) - xi(j));
  end
  y = y + d(i) * w;
  w = w .* (x - xi(i));
end
%!test
%! xi = linspace (0, 1, 11)';
%! xi(:, 2) = exp (xi);
%! x = linspace (0, 1, 101)';
%! y = intpolyvalnewt (x, xi);
%! assert (y, exp (x), 1e-12)
