function A = poisson_matrix(integerN, x0, xN)
  %
  % In the numerical solution of the poisson equation,
  %
  %   -u''(x) = f(x)
  %
  % in one dimension, subject to the boundary conditions,
  %
  %   u(x0) = 0
  %   u(xN) = 1
  %
  % over the interval [x0,xN], we need to compute a matrix A which
  % is then multiplied by the vector of u(x0), ..., u(xN). The entries
  % of A are determined from the second order finite difference formula,
  % i.e. the second order forward Euler method.
  %
  % INPUTS:
  %
  %   * ``integerN`` - An integer representing the number of
  %     subintervals we should use to approximate `u`. Must be
  %     greater than or equal to 2, since we have at least two
  %     values for u(x0) and u(xN).
  %
  %   * ``f`` - The function on the right hand side of the poisson
  %     equation.
  %
  %   * ``x0`` - The initial point.
  %
  %   * ``xN`` - The terminal point.
  %
  % OUTPUTS:
  %
  %   * ``A`` - The (N+1)x(N+1) matrix of coefficients for u(x0),
  %             ..., u(xN).

  if (integerN < 2)
    A = NA
    return
  end

  [xs,h] = partition(integerN, x0, xN);

  % We cannot evaluate u_xx at the endpoints because our
  % differentiation algorithm relies on the points directly to the left
  % and right of `x`.
  differentiable_points = xs(2:end-1);

  % These are the coefficient vectors for the u(x0) and u(xn)
  % constraints. There should be N zeros and a single 1.
  the_rest_zeros = zeros(1, integerN);
  u_x0_coeffs = cat(2, 1, the_rest_zeros);
  u_xN_coeffs = cat(2, the_rest_zeros, 1);

  % Start with the u(x0) row.
  A = u_x0_coeffs;

  for x = differentiable_points
    % Append each row obtained from the forward Euler method to A.
    u_row = forward_euler(2, xs, x);
    A = cat(1, A, u_row);
  end

  % Finally, append the last row for xN and negate the whole thing (see
  % the definition of the poisson problem).
  A = - cat(1, A, u_xN_coeffs);
end
