deltaxsq = deltax**2
deltaysq = deltay**2
DOUBLE PRECISION :: a, b, c, d, max, current
a = deltaxsq * deltaysq
b = acoeff * deltaysq
c = bcoeff * deltaxsq
d = -2*(b + c)

do m=1,5
    uitm = ufct
    max = 0
    do j=1, jmax

        do k=1, kmax

            current = (rhs(j,k)*a - b*(ufct(j-1,k)+ufct(j+1,k)) - c*(ufct(j,k-1)+ufct(j,k+1)))/d
            if(uitm(j,k)-current > max) then
                max = uitm(j,k)-current
            ufct(j,k) = current

        end do

    end do

    if(max .le. epsilon)
        exit

end do