Solving Ax=b with matrices (bigger than 2 by 2 matrices)

Question

I have been working on creating a program that solves linear systems of equations for the Jacobi and Gauss-Seidel iterative methods. (Link to the methods: https://www.cis.upenn.edu/~cis515/cis515-12-sl5.pdf) The specific problem I have with my code is that it works for 2 by 2 matrices but does not work for anything larger than that. My sample 2x2 matrix is A=[[3,1],[2,-1]] and b=[5,0]. This works just fine but the 3x3 matrix I am trying to solve is A=[[1,1,1],[0,2,5],[2,5,-1]] and b=[6,-4,27] and it does not compute the correct answer.

import numpy as np
from math import *
from matplotlib import pyplot as plt
x_list = []
x_list_1 = []

it_1 = int(input("What is the maximum amount of iterations: "))
it= it_1+1
n = int(input("Enter the matrix size: "))
print ("Enter your matrix A:")

A= np.zeros((n,n))

for i in range(n):
   A[i] = input().split(" ")

MatrixA = np.array(A)
print("A= ")
print (A)

for i in range(A.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]

print("Enter your solution b:")
b = np.zeros((n,1))

for i in range(n):
   b[i] = input().split(" ")

MatrixB = np.array(b)
print("b= ")
print (b)

for i in range(b.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(b.shape[1])]

x = np.zeros_like(b)
for it_count in range(it):
    print("Iteration", it_count, ":")
    print(x)
    x_new = np.zeros_like(x)

    for i in range(A.shape[0]):
        norm = np.linalg.norm(x)
        x_list.append(norm)
        s1 = np.dot(A[i, :i], x[:i])
        s2 = np.dot(A[i, i + 1:], x[i + 1:])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]

    if np.allclose(x, x_new, atol=1e-25, rtol=0):
        break

    x = x_new

print("Solution (Jacobi):")
print(x)

x = np.zeros_like(b)
def gauss(A, b, x, it):

    counting = range(0,it)
    lower_triangle = np.tril(A)
    upper_triangle = A - lower_triangle
    for it_count in range(it):
        print("Iteration", it_count, ":")
        print(x)
        norm = np.linalg.norm(x)
        x_list_1.append(norm)
        x = np.dot(np.linalg.inv(lower_triangle), b - np.dot(upper_triangle, x))
    return x


print("")
z = gauss(A, b, x, it)
print("Solution(Gauss-Seidel):")
print(z)

Jacobi_Magnitude_List = list(dict.fromkeys(x_list))
plt.plot(Jacobi_Magnitude_List, label="Jacobi Magnitude")
plt.plot(x_list_1, label = "Gauss-Seidel Magnitude")
plt.title("Jacobi and Gauss-Seidel Method Magnitude at each Iteration")
plt.xlabel("Iteration Number")
plt.ylabel("Magnitude")
plt.legend()
plt.show()

I am not sure what in my code is limiting it to only solving 2 by 2 matrices.