function u = solve_poisson(integerN, f)
  ##
  ## Numerically solve the poisson equation,
  ##
  ##   -u_xx(x) = f(x)
  ##
  ## in one dimension, subject to the boundary conditions,
  ##
  ##   u(0) = 0
  ##   u(1) = 1
  ##
  ## over the interval [0,1]. It is assumed that the function `f` is at
  ## least continuous on [0,1].
  ##
  ## INPUTS:
  ##
  ##   * ``integerN`` - An integer representing the number of subintervals
  ##     we should use to approximate `u`.
  ##
  ##   * ``f`` - The function on the right hand side of the poisson
  ##     equation.
  ##
  ## OUTPUTS:
  ##
  ##   * u - The function `u` that approximately solves -u_xx(x) = f(x).
  ##

  [xs,h] = partition(integerN, 0, 1);

  ## 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)

  ## Prepend and append a zero row to f(xs).
  ## Oh, and negate the whole thing (see the formula).
  F = - cat(2, 0, f(differentiable_points), 0)

  ## 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_x4_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 x4.
  A = cat(1, A, u_x4_coeffs);

  ## Solve the system. This gives us the value of u(x) at xs.
  U = A \ F';

  plot(xs,U)
end
