The sum of the second- and fourth-order Hylleraas functionals and the Newton-Raphson optimization technique have been used to optimize a manifold of first- and second-order correlation orbitals (FOCOs and SOCOs) for the spin-restricted Hartree-Fock (RHF) wave function. The correlation orbitals are linear combinations of the RHF virtual orbitals. The computational implementation of the new method is discussed and a numerical application to a simple model system is presented to illustrate the capabilities of the proposed methodology.