that it is possible to derive the solution function by a simple inversion, but this doesn't check that all eigenvalues are larger than one. It is better to use the standard algorithm with QZ decomposition