We want to solve the equation $Ax = b$. Recall that every root-finding problem has equivalent, associated fixed-point iterations. If $x$ is a solution of $Ax = b$ then $x$ satisfies
$$x = (I - A)x + b.$$and we want to find a fixed-point of
$$g(x) = (I - A)x + b.$$import numpy as np
A = np.array([[1.1, .1, .3, .4],
[.2, .9, -.1, .33],
[-.3, .4, 1., .8],
[.1, .2, .3, 1]])
b = np.array([[1],
[1],
[1],
[1]])
T = np.eye(4) - A
x = np.array([[0],
[0],
[0],
[0]])
for i in range(40):
x = np.dot(T,x) + b
np.dot(A,x) - b
array([[ -5.30764743e-09], [ -2.58408872e-09], [ -3.60924202e-09], [ -3.46897711e-09]])