The power method for solving N-particle eigenvalue equations is contracted onto the two-particle space to produce a reduced "variational" method for solving the contracted Schrodinger equation (CSE), also known as the density equation. In contrast to the methods which solve a system of approximate nonlinear equations to determine the two-particle reduced density matrix (2-RDM) nonvariationally, the contracted power method updates the 2-RDM iteratively through a "gradient" of the N-particle energy. After each power iteration we modify the 2-RDM to satisfy certain N-representability conditions through an extension of purification to correlated RDMs. The contracted power method is illustrated with a variety of molecules. Significant features of the present calculations include (i) accurate results for both first- and second-order functionals for building the 3- and the 4-RDM's from the 2-RDM's; (ii) the first molecular implementation of the Mazziotti correction within the CSE [Mazziotti, Phys. Rev. A 60, 3618 (1999)]. (iii) a spin-orbital formulation; (iv) the treatment of both core and valence orbitals as active. and; (v) a reduction of the CSE computational scaling through fast summation and the natural-orbital transformation.