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   %
   7   % INPUTS:
   8   %
   9   %   * ``f`` - The function whose divided differences we want.
  10   %
  11   %   * ``xs`` - A vector containing x-coordinates. The length of `xs`
  12   %   determines the order of the divided difference.
  13   %
  14   %
  15   % OUTPUTS:
  16   %
  17   %   * ``dd`` - The divided difference f[xs(1), xs(2),...]
  18   %
  19
  20   order = length(xs) - 1;
  21
  22   if (order < 0)
  23     % Can't do anything here. Return nothing.
  24     dd = NA;
  25   elseif (order == 0)
  26     % Our base case.
  27     dd = f(xs(1));
  28   else
  29     % Order >= 1.
  30     cs = divided_difference_coefficients(xs);
  31     dd = dot(cs, f(xs));
  32   end
  33 end