Quasi-exactly solvable models based on special functions

We suggest a systematic method of extension of quasi-exactly solvable (QES) systems. We construct finite-dimensional subspaces on the basis of special functions (hypergeometric, Airy, Bessel ones) invariant with respect to the action of differential operators of the second order with polynomial coefficients. As a example of physical applications, we show that the known two-photon Rabi Hamiltonian becomes quasi-exactly solvable at certain values of parameters when it can be expressed in terms of corresponding QES operators related to the hypergeometric function.