#!/usr/bin/octave --silent

function dd = divided_difference(f, xs)
  ## Compute divided difference of `f` at points `xs`. The argument `xs`
  ## is assumed to be a vector containing at least one element. If it
  ## contains n elements, the (n-1)st divided difference will be
  ## calculated.
  order = length(xs);
  
  if (order < 1)
    ## Can't do anything here. Return nothing.
    dd = NA;
  elseif (order == 1)
    ## Our base case.
    dd = f(xs(1))
  else
    ## Order > 1, recurse.

    ## f[x0,...,x_n-1]
    f0 = divided_difference(f, xs(1:end-1));
    ## f[x1,...,x_n]
    f1 = divided_difference(f, xs(2:end));

    # http://mathworld.wolfram.com/DividedDifference.html
    dd = (f0 - f1)/(xs(1) - xs(end))
  end
end
