A two-stage iterative scheme is proposed to handle a central problem of molecular dynamics, the computation of interior eigenvalues of large Hamiltonian matrices. The proposed method involves an initial propagation process for a time-dependent wave operator which is then inserted in an iterative process (recursive distorted wave approximation or single cycle method) to yield the exact stationary wave operator. The merits of the wave operator formalism for quasiadiabatic propagation are analyzed, and possible improvements such as the use of partial adiabatic representations and spectral filters, are outlined. The proposed algorithm is applied to the test case of two coupled oscillators with variable coupling strength, and yields accurate results even with small switching times.