Square-root solver for discrete-time Lyapunov equations
R = dlyapchol(A,B)
X = dlyapchol(A,B,E)
R = dlyapchol(A,B) computes a Cholesky factorization X = R'*R of the solution X to the Lyapunov matrix equation:
A*X*A'- X + B*B' = 0
All eigenvalues of A matrix must lie in the open unit disk for R to exist.
X = dlyapchol(A,B,E) computes a Cholesky factorization X = R'*R of X solving the Sylvester equation
A*X*A' - E*X*E' + B*B' = 0
All generalized eigenvalues of (A,E) must lie in the open unit disk for R to exist.
 Bartels, R.H. and G.W. Stewart, "Solution of the Matrix Equation AX + XB = C," Comm. of the ACM, Vol. 15, No. 9, 1972.
 Hammarling, S.J., "Numerical solution of the stable, non-negative definite Lyapunov equation," IMA J. Num. Anal., Vol. 2, pp. 303-325, 1982.
 Penzl, T., "Numerical solution of generalized Lyapunov equations," Advances in Comp. Math., Vol. 8, pp. 33-48, 1998.