solver.py

import numpy as np
def solve(A, b):
    A=np.array(A,float)
    b=np.array(b,float)
    n=len(b)
    for k in range(0,n-1):
        for i in range(k+1,n):
            if (A[i,k] != 0.0) :
                lam=A[i,k]/A[k,k]
                A[i,k]=0
                A[i,k+1:n]-=lam*A[k,k+1:n]
                b[i]-=lam*b[k]
    x=b.copy()
    for k in range(n-1,-1,-1):
        x[k]=(b[k]-np.dot(A[k,k+1:n],x[k+1:n]))/A[k,k]
    print (x)
A = [
[4, -2, 1],
[-2, 4, -2],
[1, -2, 4]
]
b = [
11,
-16,
17
]
solve(A, b)