We propose a method for solving the nuclear eigenvalue problem using a Hamiltonian of general form in a space spanned by multiphonon states built out of particle-hole TammDancoff phonons. The method consists in deriving for a given n-phonon subspace a set of equations of motion of simple structure, where all quantities are expressed in terms of energies and one-body density matrices determined in the (n − 1)-phonon subspace. Their solution, therefore, yields states built out of particle-hole operators applied on top of the (n−1)phonons basis states. Because of this peculiar structure, the states so generated form an overcomplete set. The redundancy, however, is removed completely and efficiently by a procedure based on the Choleski decomposition. The phonon composition of the basis states allows to remove naturally and maximally the spurious admixtures induced by the center of mass motion. The iteration of the equations from n = 0 up to an arbitrary number N of phonons generates an orthogonal basis spanning a space which is the direct sum of 0+1+..+N phonon subspaces. This basis decomposes the full Hamiltonian into N diagonal blocks plus off-diagonal terms coupling the n-phonons with the (n ± 1)and (n ± 2)-phonon subspaces. These off-diagonal terms are easily computed by recursive relations. The Hamiltonian can, therefore, be easily diagonalized by resorting to importance sampling techniques [1,2], which allow to truncate the multiphonon space while monitoring the accuracy of the solutions. The method can be easily upgraded so as to generate a particle-hole or quasiparticle RPA multiphonon basis. It is suitable for describing with high accuracy low-lying multiphonon spectra, whose experimental evidence is rapidly growing, as well as high-energy resonances like the recently discovered double dipole giant resonance. It can account with great accuracy for the anharmonicities associated to all collective modes. Last, but not least, it yields a highly correlated ground state which contains explicitly up to N particle N hole configurations and, therefore, applies to the cases where the quasiboson approximation breaks down. For illustrative purposes, we apply the method to 16O, whose low-lying states are dominated by complex configurations containing up to 4 particle 4 hole states. Such a calculation probes the extreme accuracy of the method, on the one hand, and, on the other hand, provides a complete and exhaustive information on the microscopic structure of low and high energy states in 16O.