function coefficients = central_difference(xs, x)
  %
  % The first order central difference at x1 is,
  %
  %   f'(x1) = (f(x2) - f(x0))/2
  %
  % where the index x1 is of course arbitrary but x2, x0 are adjacent
  % to x1. The coefficients we seek are the coefficients of f(xj) for
  % j = 1,...,N-2, where N is the length of ``xs``. We omit the first
  % and last coefficient because at x0 and xN, the previous/next
  % value is not available.
  %
  % This should probably take an 'order' parameter as well; see
  % forward_euler().
  %
  % INPUT:
  %
  %   * ``xs`` - The vector of x-coordinates.
  %
  %   * ``x`` - The point `x` at which you'd like to evaluate the
  %   derivative of the specified `integer_order`. This should be an
  %   element of `xs`.
  %
  % OUTPUT:
  %
  %   * ``coefficients`` - The vector of coefficients, in order, of
  %   f(x0), f(x1), ..., f(xn).
  %

  if (length(xs) < 3)
    % We need at least one point other than the first and last.
    coefficients = NA;
    return;
  end

  x_idx = find(xs == x);

  if (x_idx == 1 || x_idx == length(xs))
    % You asked for the difference at the first or last element, which
    % we can't do.
    coefficients = NA;
    return;
  end

  % Start with a vector of zeros.
  coefficients = zeros(1, length(xs));

  % And fill in the two values that we know.
  coefficients(x_idx - 1) = -1/2;
  coefficients(x_idx + 1) = 1/2;
end
