l = 40;
n = l*l;

r = zeros(n,1);
r(idx(l/2-l/5,l/2,l)) = -1;
r(idx(l/2+l/5,l/2,l)) = +1;

A = eye(n);

for i = 2:l-1
  for j = 2:l-1
    k = idx(i,j,l);
    A(k,k) = -4;
    A(k,k+1) = 1;
    A(k,k-1) = 1;
    A(k,k+l) = 1;
    A(k,k-l) = 1;
  end
end

u = A\r;

surfc(reshape(u,l,l));

