In neutron chain systems with material symmetries, various k-eigenvalues of the neutron balance equation beyond the dominant one may be degenerated. As shown in a companion paper, the power iteration method can be used to compute higher eigenfunctions in symmetric systems, provided that the global problem is partitioned into symmetry class–related lower-sized problems with appropriate boundary conditions. Those boundary conditions have been implemented in the diffusion solver of the ERANOS code system in rectangular geometry and within the framework of a discontinuous Galerkin spatial approximation of the multigroup discrete ordinates transport equation in the SNATCH solver. Numerical results in homogeneous geometry are provided for verification purposes, as well as the first eigenfunctions of the Takeda benchmarks. Finally, the transport effect on the first flux harmonics for an industrial-sized reactor ZPPR-18 is discussed.