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

ts = clock();

l = 20;
n = l*l;

r = zeros(n,1);
r(idx(l/2-3,l/2,l)) = -1;
r(idx(l/2+3,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;
  endfor
endfor

u = A\r;

etime(clock(),ts)

surf(reshape(u,l,l));

