It is need to solve the problem using the Jacobi method: d^T/dx^2 + d^2T/dy^2 = 0 with conditions at the boundaries of the square: T(0, y) = 1, T(x, 1) = 1, dT/dx(1, y) = 1 - y, dT/dy(x, 0) = x.

Conditions for stopping iterations: max|T_new - T| < d_t and max|R_{i, j}| < d_R, where R is the residual matrix.

I approximate the boundary conditions on derivatives by 2-order finite-difference schemes. By the condition of the problem, the step along the x and y axes is h = 1/20. My implementation of the numerical method results in 645 iterations for convergence in the case of the Jacobi method and 999 in the case of Gauss-Seidel, which is most likely incorrect. I think that the problem is in the Jacobi method. Please help me find out why this can happen. I enclose my code for the Jacobi method.

h = 1/20
grid = np.arange(0, 1+h, h)
N = len(grid)
 
def jacobi_method(M):
    T  = M.copy() 
    T[0, :] = 1
    T[:, -1] = 1
    T_prev = T.copy()
    for i in range(1, N - 1):
        for j in range(1, N - 1):
            T[i, j] = ( T_prev[i+1, j] + T_prev[i-1, j] + T_prev[i, j+1] + T_prev[i, j-1]) /4
    
    T[-1, :] = (4*T[-2, : ] - T[-3, :] + 2*h*(1- grid))/3
    T[:,0] = (4*T[:, 1] - T[:, 2] - 2*h*grid)/3
 
    return T
 
def residual(T):
    residual = np.zeros((len(T)-2, len(T)-2))
    for i in range(1, len(T) - 1):
        for j in range(1, len(T[i]) - 1):
            residual[i-1, j-1] = ((T[i+1, j] - 2*T[i, j] + T[i-1, j]) / (h**2) +
                                  (T[i, j+1] - 2*T[i, j] + T[i, j-1]) / (h**2))
 
    return residual
    
T1 = np.ones((N, N))
dT = 1
deltaT = 10e-7
R = 1
deltaR = 10e-5
iterations = []
dTs = []
Rs = []
i = 0
while dT >= deltaT or R >= deltaR :
    T_n = jacobi_method(T1)
    dT = np.max(abs(T_n - T1))
    dTs.append(dT)
    R = np.max(abs(residual(T_n)))
    Rs.append(R)
    T1 = T_n.copy()
    i += 1
    iterations.append(i)