compute the residual at k-th step using the new solution
for j = 2:ny for i = 2:nx
d2udx2 = (un[i+1,j] - 2*un[i,j] + un[i-1,j])/(dx^2)
d2udy2 = (un[i,j+1] - 2*un[i,j] + un[i,j-1])/(dy^2)
r[i,j] = f[i,j] - d2udx2 - d2udy2
end~end