def QRsolve(self, b):

        Q, R = self.as_mutable().QRdecomposition()
        y = Q.T*b

        # back substitution to solve R*x = y:
        # We build up the result "backwards" in the vector 'x' and reverse it
        # only in the end.
        x = []
        n = R.rows
        for j in range(n - 1, -1, -1):
            tmp = y[j, :]
            for k in range(j + 1, n):
                tmp -= R[j, k]*x[n - 1 - k]
            x.append(tmp / R[j, j])
        return self._new([row._mat for row in reversed(x)])
Solve the linear system 'Ax = b'.
dataset/ETHPy150Open sympy/sympy/sympy/matrices/matrices.py/MatrixBase.QRsolve