function A = laplacian2d(n)
% LAPLACIAN2D - generate sparse 2D discrete Laplacian matrix of dimension
%               n^2 x n^2, on a square domain.
    r = zeros(1, n);
    r(1:2) = [2, -1];   
    T = toeplitz(r);
    E = speye(n);                % nxn identity matrix, stored sparsely
    A = kron(T, E) + kron(E, T); % Kronecker sum, stored sparsely
end
