Проблема с решением системы линейных алгебраических уравнений методом нижней релаксации

Всем добрый день. Нужна ваша помощь. По заданию нужно решить СЛАУ методом нижней релаксации. Но мой код работает только с матрицами до 3х3, а начиная с 4х4 работает только с очень маленькими значениями. И даже если задать маленькое количество итераций для матрицы 4х4, то в результате выдает значения, которые даже близко не похожи на истинные. Помогите пожалуйста...

import numpy as np

def relax_solver(A, b, omega, initial_guess, convergence_criteria, stop, itera):
    """
    Solve the system of linear equations using the lower relaxation method.
    Inputs:
        A: nxn numpy matrix
        b: n dimensional numpy vector
        omega: relaxation factor
        initial_guess: an initial solution guess for the solver to start with
        stop: maximum number of iterations.
        itera: a flag that determines whether information about each iteration should be printed
    Returns:
        None
    """
    phi = initial_guess[:]
    residual = np.linalg.norm(np.matmul(A, phi) - b) 
    iteration = 0
    while residual > convergence_criteria and stop > iteration:
        for i in range(A.shape[0]):
            sigma = 0
            for j in range(A.shape[1]):
                if j != i:
                    sigma += A[i][j] * phi[j]
            phi[i] = (1 - omega) * phi[i] + (omega / A[i][i]) * (b[i] - sigma)
        residual = np.linalg.norm(np.matmul(A, phi) - b)
        iteration += 1
        if itera == True:
            print(f'Iteration {iteration}, Solution: {", ".join(f"{val:.15f}" for val in phi)}')
            
    print("Solution:", ", ".join(f"{val:.15f}" for val in phi))
    print(f"Number of iterations: {iteration}")

Пробовал с такими входными данными:

A = np.array([[5, 8, -3, 21],
              [-3, 4, 0.5, -8.7],
              [-3.6, 0.34, 3, -76],
              [-5.7, 2.26, -5.6, 9]])

B = np.array([-6, 4, 45, -15], dtype=float)