ts = clock();

n = 20;
u = zeros(n);
r = zeros(n);

r(n/2-3,n/2) = -1;
r(n/2+3,n/2) = +1;

for t = 0:100
  for i = 2:n-1
    for j = 2:n-1
      u(i,j) = 0.25*(u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - r(i,j));
    endfor
  endfor
endfor

etime(clock(),ts)

surf(u);
