Working within an exactly solvable 3 level model, we discuss am extension of the Random Phase Approximation (RPA) based on a boson formalism. A boson Hamiltonian is defined via a mapping procedure and its expansion truncated at four-boson terms. RPA-type equations are then constructed and solved iteratively. The new solutions gain in stability with respect to the RPA ones. We perform diagonalizations of the boson Hamiltonian in spaces containing up to four-phonon components. Approximate spectra exhibit an improved quality with increasing the size of these multiphonon spaces. Special attention is addressed to the problem of the anharmonicity of the spectrum.