We present a novel and simple algorithm in the variation after projection (VAP) approach for the non-yrast nuclear states. It is for the first time that the yrast state and non-yrast states can be varied on the same footing. The orthogonality among the calculated states is automatically fulfilled by solving the Hill-Wheeler equation. This avoids the complexity of the frequently used Gram-Schmidt orthogonalization, as adopted by the excited VAMPIR method. Thanks to the Cauchys interlacing theorem in the matrix theory, the sum of the calculated lowest projected energies with the same quantum numbers can be safely minimized. Once such minimization is converged, all the calculated energies and the corresponding states can be obtained, simultaneously. The present VAP calculations are performed with time-odd Hartree-Fock Slater determinants. It is shown that the calculated VAP energies (both yrast and non-yrast) are very close to the corresponding ones from the full shell model calculations. It looks the present algorithm is not limited to the VAP, but should be universal, i.e., one can do the variation with different forms of the many-body wavefunctions to calculate the excited states in different quantum many-body systems.