From 6de25e4c4273789fd2f98e309db4ef3f0902b185 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 15 Oct 2012 15:00:55 -0400
Subject: [PATCH] Factor out the poisson_matrix function.

---
 poisson_matrix.m | 65 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 65 insertions(+)
 create mode 100644 poisson_matrix.m

diff --git a/poisson_matrix.m b/poisson_matrix.m
new file mode 100644
index 0000000..696185f
--- /dev/null
+++ b/poisson_matrix.m
@@ -0,0 +1,65 @@
+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.
+  A = cat(1, A, u_xN_coeffs);
+end
-- 
2.43.2