function coefficients = forward_euler(integer_order, xs, x)
  %
  % Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
  % Take for example a first order approximation, with,
  %
  %   xs = [x0,x1,x2,x3,x4]
  %
  %   f'(x1) ~= [f(x2)-f(x1)]/(x2-x1)
  %
  % This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids the
  % solution of linear systems.
  %
  %
  % INPUTS:
  %
  %   * ``integer_order`` - The order of the derivative which we're
  %     approximating.
  %
  %   * ``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`.
  %
  %
  % OUTPUTS:
  %
  %   * ``coefficients`` - The vector of coefficients, in order, of
  %     f(x0), f(x1), ..., f(xn).
  %

  if (integer_order < 0)
    % You have made a grave mistake.
    coefficients = NA;
    return;
  end

  if (integer_order == 0)
    coefficients = x;
    return;
  end

  if (length(xs) < 2)
    % You can't approximate a derivative of order greater than zero
    % with zero or one points!
    coefficients = NA
    return;
  end

  if (even(integer_order))
    offset_b = integer_order / 2;
    offset_f = offset_b;
  else
    % When the order is odd, we need one more "forward" point than we
    % do "backward" points.
    offset_b = (integer_order - 1) / 2;
    offset_f = offset_b + 1;
  end

  % Zero out the coefficients for terms that won't appear. We compute
  % where `x` is, and we just computed how far back/forward we need to
  % look from `x`, so we just need to make the rest zeros.
  x_idx = find(xs == x);
  first_nonzero_idx = x_idx - offset_b;
  last_nonzero_idx = x_idx + offset_f;
  leading_zero_count = first_nonzero_idx - 1;
  leading_zeros = zeros(1, leading_zero_count);
  trailing_zero_count = length(xs) - last_nonzero_idx;
  trailing_zeros = zeros(1, trailing_zero_count);

  targets = xs(first_nonzero_idx : last_nonzero_idx);

  % The multiplier comes from the Taylor expansion.
  multiplier = factorial(integer_order);
  cs = divided_difference_coefficients(targets) * multiplier;

  coefficients = horzcat(leading_zeros, cs, trailing_zeros);
end
