--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/neumann.m	Thu May 07 18:36:24 2015 +0200
@@ -0,0 +1,24 @@
+function [A, T] = neumann(n)
+%NEUMANN  Singular matrix from the discrete Neumann problem (sparse).
+%         NEUMANN(N) is the singular, row diagonally dominant matrix resulting
+%         from discretizing the Neumann problem with the usual five point
+%         operator on a regular mesh.
+%         It has a one-dimensional null space with null vector ONES(N,1).
+%         The dimension N should be a perfect square, or else a 2-vector,
+%         in which case the dimension of the matrix is N(1)*N(2).
+
+%         Reference:
+%         R.J. Plemmons, Regular splittings and the discrete Neumann
+%         problem, Numer. Math., 25 (1976), pp. 153-161.
+
+if max(size(n)) == 1
+   m = sqrt(n);
+   if m^2 ~= n, error('N must be a perfect square.'), end
+   n(1) = m; n(2) = m;
+end
+
+T = tridiag(n(1), -1, 2, -1);
+T(1,2) = -2;
+T(n(1),n(1)-1) = -2;
+
+A = kron(T, eye(n(2))) + kron(eye(n(2)), T);