def LUsolve(self, rhs, iszerofunc=_iszero):
        if rhs.rows != self.rows:
            raise ShapeError("`self` and `rhs` must have the same number of rows.")

        A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero)
        n = self.rows
        b = rhs.permuteFwd(perm).as_mutable()
        # forward substitution, all diag entries are scaled to 1
        for i in range(n):
            for j in range(i):
                scale = A[i, j]
                b.zip_row_op(i, j, lambda x, y: x - y*scale)
        # backward substitution
        for i in range(n - 1, -1, -1):
            for j in range(i + 1, n):
                scale = A[i, j]
                b.zip_row_op(i, j, lambda x, y: x - y*scale)
            scale = A[i, i]
            b.row_op(i, lambda x, _: x/scale)
        return rhs.__class__(b)
Solve the linear system Ax = rhs for x where A = self.
dataset/ETHPy150Open sympy/sympy/sympy/matrices/matrices.py/MatrixBase.LUsolve