function S = advection_matrix(integerN, x0, xN)
  ##
  ## The numerical solution of the advection-diffusion equation,
  ##
  ##   -d*u''(x) + v*u'(x) + r*u = f(x)
  ##
  ## in one dimension, subject to the boundary conditions,
  ##
  ##   u(x0)  = u(xN)
  ##
  ##   u'(x0) = u'(xN)
  ##
  ## over the interval [x0,xN] gives rise to a linear system:
  ##
  ##   AU = h^2*F
  ##
  ## where h = 1/n, and A is given by,
  ##
  ##   A = d*K + v*h*S + r*h^2*I.
  ##
  ## We will call the matrix S the "advection matrix," and it will be
  ## understood that the first row (corresponding to j=0) is to be
  ## omitted; since we have assumed that when j=0, u(xj) = u(x0) =
  ## u(xN) and likewise for u'. ignored (i.e, added later).
  ##
  ## 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:
  ##
  ##   * ``S`` - The NxN matrix of coefficients for the vector [u(x1),
  ##   ..., u(xN)].
  ##
  ## EXAMPLES:
  ##
  ## For integerN=4, x0=0, and x1=1, we will have four subintervals:
  ##
  ##   [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]
  ##
  ## The first row of the matrix 'S' should compute the "derivative"
  ## at x1=0.25. By the finite difference formula, this is,
  ##
  ##   u'(x1) = (u(x2) - u(x0))/2
  ##
  ##          = (u(x2) - u(x4))/2
  ##
  ## Therefore, the first row of 'S' should look like,
  ##
  ##   2*S1 = [0, 1, 0, -1]
  ##
  ## and of course we would have F1 = [0] on the right-hand side.
  ## Likewise, the last row of S should correspond to,
  ##
  ##   u'(x4) = (u(x5) - u(x3))/2
  ##
  ##          = (u(x1) - u(x3))/2
  ##
  ## So the last row of S will be,
  ##
  ##   2*S4 = [1, 0, -1, 0]
  ##
  ## Each row 'i' in between will have [-1, 0, 1] beginning at column
  ## (i-1). So finally,
  ##
  ##   2*S = [0, 1, 0, -1]
  ##	     [-1, 0, 1, 0]
  ##         [0, -1, 0, 1]
  ##         [1, 0, -1, 0]

  if (integerN < 2)
    S = 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`. Since we're starting at j=1 anyway, we cut
  ## off two from the beginning.
  differentiable_points = xs(3: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 - 3);
  u_x0_coeffs = cat(2, the_rest_zeros, [0.5, 0, -0.5]);
  u_xN_coeffs = cat(2, [0.5, 0, -0.5], the_rest_zeros);

  ## Start with the u(x0) row.
  S = u_x0_coeffs;

  for x = differentiable_points
    ## Append each row obtained from the forward Euler method to S.
    ## Chop off x0 first.
    u_row = central_difference(xs(2:end), x);
    S = cat(1, S, u_row);
  end

  ## Finally, append the last row for xN.
  S = cat(1, S, u_xN_coeffs);
end
