% interpolate.m
%
% example usages:
%
%   > x = [0 1 2]
%   > y = [0 1 4]
%   > t = linspace(-1, 3, 100)
%   > s = interpolate(x, y, t)
%   > plot(x, y, 'ro', t, s, 'g-')
%
% or
%
%   > x = linspace(-1, 1, 10) <-- try gradually increasing this 10 and see what happens!
%   > f = inline('1 ./ (1 + 12*x.^2)');
%   > t = linspace(-1, 1, 1000)
%   > s = interpolate(x, f(x), t);
%   > plot(x, f(x), 'ro', t, s, 'g-');
%

function s = interpolate(x, y, t)

  n = length(x);

  % first compute divided difference table D = [ D(i, j) ] so that
  %   D(i, j) = f[x(i), x(i+1), ..., x(i+j-1)] 
  % and then compute
  %   s = f[x(1)] + f[x(1), x(2)] * (t - x(1)) +
  %                         ... + f[x(1), ..., x(n)] * (t-x(1)) * ... * (t-x(n-1))
  % using nested multiplication.

  D = zeros(n, n);

  for c = 1:n
    for r = 1:n-c+1
      if c == 1
	D(r, c) = y(r);
      else
	D(r, c) = (D(r+1, c-1) - D(r, c-1)) / (x(r+c-1) - x(r));
      end
    end
  end
  
  s = D(1, n) + 0 * t;
  for c = n-1:-1:1
    s = D(1, c) + ((t - x(c)) .* s);
  end

end