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.
##
## INPUTS:
##
## * ``f`` - The function whose divided differences we want.
##
## * ``xs`` - A vector containing x-coordinates. The length of `xs`
## determines the order of the divided difference.
##
##
## OUTPUTS:
##
## * ``dd`` - The divided difference f[xs(1), xs(2),...]
##
order = length(xs) - 1;
if (order < 0)
## Can't do anything here. Return nothing.
dd = NA;
elseif (order == 0)
## Our base case.
dd = f(xs(1));
else
## Order >= 1.
cs = divided_difference_coefficients(xs);
dd = dot(cs, f(xs));
end
end