function d = diffdiv2 (xi)
%
% function d = diffdiv2 (xi)
%
% Calcolo delle differenze divise dei punti [xi(:,1),xi(:,2)]. I nodi
% di interpolazione possono avere molteplicita` maggiore di uno. In
% tal caso, nodi uguali devono essere consecutivi xi(m(i-1)+1:m(i),1)
% e i corrispondenti valori di xi(m(i-1)+1:m(i),2) devono contenere i
% valori della funzione da interpolare e delle successive derivate nel
% punto xi(m(i-1)).
yi = xi(:, 2);
xi = xi(:, 1);
n = length (xi) - 1;
d = zeros (n + 1,1);
t(1) = yi(1);
d(1) = t(1);
for i = 1:n + 1
  k = 0;
  t(i) = yi(i);
  for j = i - 1:-1:1
    if (xi(j) ~= xi(i))
      t(j) = (t(j + 1) - t(j)) / (xi(i) - xi(j));
    else
      k = k + 1;
      t(j) = yi(i) / factorial (k);
    end
  end
  d(i) = t(1);
end
%!demo
%! n = 5;
%! xi = linspace (0, 2 * pi, n + 1)';
%! xi(:, 2) = sin (xi(:, 1));
%! x = linspace (0, 2 * pi, 101);
%! di = diffdiv2 (xi);
%! y = zeros (size (x));
%! w = ones (size (x));
%! for i = 1:n + 1
%!   y = y + di(i) * w;
%!   w = w .* (x - xi(i, 1));
%! end
%! plot (x, y, xi(:, 1), xi(:, 2), '*')
%! hold on
%! xfin = linspace (6, 7, 51);
%! plot (xfin, (xfin - 2 * pi))
%!demo
%! n = 5;
%! xi = linspace (0, 2 * pi, n)';
%! xi = [xi; 2 * pi];
%! xi(1:n, 2) = sin (xi(1:n, 1));
%! xi(n + 1, 2) = cos (xi(n + 1, 1));
%! x = linspace (0, 2 * pi, 101);
%! di = diffdiv2 (xi);
%! y = zeros (size (x));
%! w = ones (size (x));
%! for i = 1:n + 1
%!   y = y + di(i) * w;
%!   w = w .* (x - xi(i, 1));
%! end
%! plot (x, y, xi(1:n, 1), xi(1:n, 2), '*')
%! hold on
%! xfin = linspace (6, 7, 51);
%! plot (xfin, (xfin - 2 * pi))
