
// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);

n = 20;
u = zeros(n,n);
r = zeros(n,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));
    end;
  end;
end;

surf(u);

