We present a new approach based on the Davidson algorithm which provides eigenvalues and eigenvectors of selected highly excited (ro)vibrational states of polyatomic molecules. The key ingredient is a prediagonalization–perturbation scheme applied to a subspace of a curvilinear normal modes basis set (including diagonal anharmonicities). The efficiency of this method is illustrated by computing all vibrational states of the H<sub>2sub>CO molecule, up to 9500cm<sup>−1sup> of internal excitation. Convergence of the levels can be assessed during the iteration process by looking at the residual ||(H−E<sub>αsub>)|Ψ<sub>αsub>mg align=center border=0 SRC=/images/glyphs/BEA.GIF>||.