Matrix eqn

Added solve_discrete_lyapunov and solve_discrete_riccati.
We definitely want this `A0 = (I - gamma .* G0) * A - B * (R_hat \ C)` (added parens around `I - gamma.*G0`