Often in reactor dynamics, higher eigenfunctions of the multigroup diffusion equation must be determined. An algorithm to calculate higher eigenfunctions (modes) of the λ-eigenvalue problem corresponding to the steady-state two-group neutron diffusion equation is presented. The method is based on a special type of subspace iteration for large sparse nonsymmetric eigenvalue problems. Having been tested using an International Atomic Energy Agency benchmark problem and also applied to a VVER-1000pressurized water reactor assembly, the algorithm was found to work very effectively and reliably. In its application, the algorithm presented is not restricted to the λ-eigenvalue problem only but is also generally applicable to large sparse nonsymmetric eigenvalue problems even with multiple and complex eigenvalues.