Решение методом Гаусса

def solve_lu(L, U, ans):
    n = L.shape[0]
    y = np.matrix(np.zeros([n, 1]))
    for i in range(n):
        y[i, 0] = ans[i, 0] - L[i, :i] * y[:i]

    x = np.matrix(np.zeros([n, 1]))
    for i in range(1, n + 1):
        x[-i, 0] = (y[-i] - U[-i, -i:] * x[-i:, 0]) / U[-i, -i]
    return x