calculate numerical solution at k-th step using Gauss-Seidel iterative method
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
un[i,j] = un[i,j] + r[i,j]/den # update the solution at every step
end~end