function [ A, B ] = compactlaplacian( n )
%COMPACTLAPLACIAN Fourth-order nine-point stencil Laplacian
% Mehrstellen Verfahren of Collatz (1950)
% A Del U = B U/h^2 to O(h^4)
% Periodic boundary conditions on a square equally-spaced grid
% Row-ordering of the grid
% Both A and B are block circulant with circulant blocks
% Both A and B are symmetric
% Condition number of A is bounded by 3, independent of n

E = speye(n);
C = sparse(2:n,1:n-1,1,n,n);
C(n,1)=1;
C = C + C'; % Sparse symmetric circulant
A = kron(C,E/8) + kron(E,C/8) + speye(n^2);

D = -5*E + C;
F = E + C/4;
B = kron( E, D ) + kron(C,F);

end