We study a large class of models with an arbitrary (finite) number of degrees of freedom, described by Hamiltonians which are polynomial in bosonic creation and annihilation operators, and including as particular cases n-th harmonic generation and photon cascades. For each model, we construct a complete set of commuting integrals of motion of the Hamiltonian, fully characterize the common eigenspaces of the integrals of motion, and show that the action of the Hamiltonian in these common eigenspaces can be represented by a quasi-exactly solvable reduced Hamiltonian, whose expression in terms of the usual generators of sl 2 is computed explicitly.