Refactor matrix equations (Ricatti, Lyapunov, Sylvester)

Currently these equations live in different modules and some of them are not easily discoverable.

For all of them there are also several possible implementations:

• in scipy (calling linpack)
• in quantecon (i.e. doublej)
• in Slycot, which is based on the last opensource version of SLICOT

Slycot, used to be complicated to build, but this has changed, and it is now pip-installable. It is supposed to be efficient and AFAIK is the only option for Generalized Sylvester equations.
The python-control library has a mateqn.py file (https://github.com/python-control/python-control/blob/master/src/mateqn.py) that is also quite complete and that we could partly reuse or depend on. It seems to do extensive error handling and input checks.

My point is that we should deal with all of these matrix equations for optimal control together in a systematic way. For instance we could

• move them to a common module (like matrix_equations.py)
• provide wrapper functions (solve_lyapunov, solve_generalized_sylvester, solve_riccati, etc...) that can switch between different implementations
• investigate which implementation is the best.

Start by locating the existing modules containing the Riccati, Lyapunov, and Sylvester equations, then read python-control's mateqn.py and compare the scipy, quantecon, and Slycot options described here. The work is done when the equations have a systematic common organization, wrapper functions cover the required variants, and the implementation choice is documented.