Some further ideas

Method of lines (MOL)

The reference given in the Hujanen-Kalenius documentation could be further investigated. Codes of the authors could probably be obtained. I asked their addresses but didn't yet get a reply, perhaps for instance the H-K group could try to find the codes and work on this a bit more.) See here It also revealed these Matlab tools which are certailnly worth looking into:

Matlab tools for Laplace/Poisson and gridding

>> help numgrid        

 NUMGRID Number the grid points in a two dimensional region.
    G = NUMGRID('R',n) numbers the points on an n-by-n grid in
    the subregion of -1<=x<=1 and -1<=y<=1 determined by 'R'.
    SPY(NUMGRID('R',n)) plots the points.
    DELSQ(NUMGRID('R',n)) generates the 5-point discrete Laplacian.
    The regions currently available are:
       'S' - the entire square.
       'L' - the L-shaped domain made from 3/4 of the entire square.
       'C' - like the 'L', but with a quarter circle in the 4-th square.
       'D' - the unit disc.
       'A' - an annulus.
       'H' - a heart-shaped cardioid.
       'B' - the exterior of a "Butterfly".
       'N' - a nested dissection ordering of the square.
    To add other regions, edit toolbox/matlab/demos/numgrid.m.
 
    See also DELSQ, DELSQSHOW, DELSQDEMO.

>> spy(numgrid('L',20))
>> spy(numgrid('A',20))
>> spy(numgrid('B',20))

>> help delsq

 DELSQ  Construct five-point finite difference Laplacian.
    delsq(G) is the sparse form of the two-dimensional,
    5-point discrete negative Laplacian on the grid G.
    The grid G can be generated by NUMGRID or NESTED.
 
    See also NUMGRID, DEL2, DELSQDEMO.

Moler's version of "lapmatr"

This is how we started. Here's a more advanced version using directly the sparse matrix indexing, and tools that also allow different geometries.

function D = delsq(G)
%DELSQ  Construct five-point finite difference Laplacian.
%   delsq(G) is the sparse form of the two-dimensional,
%   5-point discrete negative Laplacian on the grid G.
%   The grid G can be generated by NUMGRID or NESTED.
%
%   See also NUMGRID, DEL2, DELSQDEMO.

%   C. Moler, 7-16-91.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.3 $  $Date: 1997/11/21 23:25:31 $

[m,n] = size(G);

% Indices of interior points
p = find(G);

% Connect interior points to themselves with 4's.
i = G(p);
j = G(p);
s = 4*ones(size(p));

% for k = north, east, south, west
for k = [-1 m 1 -m]
   % Possible neighbors in k-th direction
   Q = G(p+k);
   % Index of points with interior neighbors
   q = find(Q);
   % Connect interior points to neighbors with -1's.
   i = [i; G(p(q))];
   j = [j; Q(q)];
   s = [s; -ones(length(q),1)];
end
D = sparse(i,j,s);