divided_difference.m

   1 function dd = divided_difference(f, xs)
   2   ## Compute divided difference of `f` at points `xs`. The argument `xs`
   3   ## is assumed to be a vector containing at least one element. If it
   4   ## contains n elements, the (n-1)st divided difference will be
   5   ## calculated.
   6   order = length(xs);
   7
   8   if (order < 1)
   9     ## Can't do anything here. Return nothing.
  10     dd = NA;
  11   elseif (order == 1)
  12     ## Our base case.
  13     dd = f(xs(1))
  14   else
  15     ## Order > 1, recurse.
  16
  17     ## f[x0,...,x_n-1]
  18     f0 = divided_difference(f, xs(1:end-1));
  19     ## f[x1,...,x_n]
  20     f1 = divided_difference(f, xs(2:end));
  21
  22     # http://mathworld.wolfram.com/DividedDifference.html
  23     dd = (f0 - f1)/(xs(1) - xs(end))
  24   end
  25 end