Для решения системы линейных уравнений типа Ax = b в numpy имеется функция linalg.solve(). Они принимает на вход квадратную матрицу линейно-независимых строк A, а также вектор значений b, и рассчитывает вектор корней x. Для проверки результата в данном примере использована функция allclose(), она сравнивает два массива поэлементно и возвращает True, если разница не превышает заданной величины (по-умолчанию относительная разность не более 1E-5, абсолютная - 1E-8).
Если найти решение не удалось, можно попробовать функцию linalg.lstsq(), которая ищет такой вектор x, для которого евклидова норма разности между Ax и b будет минимальна.


