Matrix eqn

Added solve_discrete_lyapunov and solve_discrete_riccati.
Do we need to make sure the matrix multiply happens before the scalar one? i.e. write `H0 = gamma .* (A' * A0) - Q_tilde`