function y = idx(i,j,nj)
  y = (i-1)*nj + (j-1) + 1;
endfunction

l = 20;
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,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;

surf(matrix(u,l,l));

