A method for solving Schrödinger’s equation that combines the Rayleigh quotient iteration (RQI) with the conjugate residual (CR) method is described. It is assumed that the eigenfunctions of the Hamiltonian operator H ˆ are expressed as linear combinations of plane waves, which allows the efficient computation of the product H ⋅ x via fast Fourier transform, where H is the representation of H ˆ in the plane-wave basis and x an arbitrary vector. The algorithm carries out two kinds of iterations: outer-loop RQI iterations (very few) and inner-loop CR iterations (many more). A very small number of orthogonalizations of the set of eigenvectors is required by the proposed method, just one per RQI iteration, which is one of its key advantages because, for large enough systems, orthogonalization is the costliest component of the calculation. The procedure is much faster than matrix diagonalization, and its cost with respect to the size of the plane-wave basis ( N P W ) has an effective scaling ranging from linear O ( N P W ) to quadratic O ( N P W 2 ) .