How to use the quantecon.matrix_eqn.solve_discrete_lyapunov function in quantecon

np.linalg.matrix_power(self.A0, 5)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])

        # === Gross rates of return on 1-period risk-free bonds === #
        self.R1_Gross = 1 / self.R1_Price

        # === Net rates of return on J-period risk-free bonds === #
        # === See p.148: log of gross rate of return, divided by j === #
        self.R1_Net = np.log(1 / self.R1_Price) / 1
        self.R2_Net = np.log(1 / self.R2_Price) / 2
        self.R5_Net = np.log(1 / self.R5_Price) / 5

        # === Value of asset whose payout vector is Pay*xt === #
        # See p.145: Equation (7.11.1)
        if isinstance(Pay, np.ndarray) == True:
            self.Za = Pay.T.dot(self.Mc)
            self.Q = solve_discrete_lyapunov(
                self.A0.T * self.beta**0.5, self.Za)
            self.q = self.beta / (1 - self.beta) * \
                np.trace(self.C.T.dot(self.Q).dot(self.C))
            self.Pay_Price = np.empty((ts_length + 1, 1))
            self.Pay_Gross = np.empty((ts_length + 1, 1))
            self.Pay_Gross[0, 0] = np.nan
            for i in range(ts_length + 1):
                self.Pay_Price[i, 0] = (xp[:, i].T.dot(self.Q).dot(
                    xp[:, i]) + self.q) / e1.dot(self.Mc).dot(xp[:, i])
            for i in range(ts_length):
                self.Pay_Gross[i + 1, 0] = self.Pay_Price[i + 1,
                                                          0] / (self.Pay_Price[i, 0] - Pay.dot(xp[:, i]))
        return

P, F, d = lq.stationary_values()

    # Need y_0 to compute government tax revenue.
    P21 = P[3, :3]
    P22 = P[3, 3]
    z0 = np.array([1, Q0, τ0]).reshape(-1, 1)
    u0 = -P22**(-1) * P21 @ z0
    y0 = np.vstack([z0, u0])

    # Define A_F and S matricies
    AF = A - B @ F
    S = np.array([0, 1, 0, 0]).reshape(-1, 1) @ np.array([[0, 0, 1, 0]])

    # Solves equation (25)
    temp = β * AF.T @ S @ AF
    Ω = solve_discrete_lyapunov(np.sqrt(β) * AF.T, temp)
    T0 = y0.T @ Ω @ y0

    return T0, A, B, F, P

P, F, d = lq.stationary_values()

    # Need y_0 to compute government tax revenue.
    P21 = P[3, :3]
    P22 = P[3, 3]
    z0 = np.array([1, Q0, tau0]).reshape(-1, 1)
    u0 = -P22**(-1) * P21.dot(z0)
    y0 = np.vstack([z0, u0])

    # Define A_F and S matricies
    AF = A - B.dot(F)
    S = np.array([0, 1, 0, 0]).reshape(-1, 1).dot(np.array([[0, 0, 1, 0]]))

    # Solves equation (25)
    temp = beta * AF.T.dot(S).dot(AF)
    Omega = solve_discrete_lyapunov(np.sqrt(beta) * AF.T, temp)
    T0 = y0.T.dot(Omega).dot(y0)

    return T0, A, B, F, P

unity
    B : array_like(float, ndim=2)
        An n x n matrix as described above.  We assume in order for
        convergence that the eigenvalues of :math:`A` have moduli bounded by
        unity
    max_it : scalar(int), optional(default=50)
        The maximum number of iterations

    Returns
    ========
    gamma1: array_like(float, ndim=2)
        Represents the value :math:`V`

    """

    gamma1 = solve_discrete_lyapunov(A, B, max_it)

    return gamma1

"""
        # == Simplify names == #
        Q, R, A, B, C = self.Q, self.R, self.A, self.B, self.C
        beta, theta = self.beta, self.theta

        # == Solve for policies and costs using agent 2's problem == #
        K_F, P_F = self.F_to_K(F)
        I = np.identity(self.j)
        H = inv(I - C.T.dot(P_F.dot(C)) / theta)
        d_F = log(det(H))

        # == Compute O_F and o_F == #
        sig = -1.0 / theta
        AO = sqrt(beta) * (A - dot(B, F) + dot(C, K_F))
        O_F = solve_discrete_lyapunov(AO.T, beta * dot(K_F.T, K_F))
        ho = (trace(H - 1) - d_F) / 2.0
        tr = trace(dot(O_F, C.dot(H.dot(C.T))))
        o_F = (ho + beta * tr) / (1 - beta)

        return K_F, P_F, d_F, O_F, o_F