The spherical harmonics or PN method is intended to approximate the neutron angular flux by a linear combination of spherical harmonics of degree at most . In this work, the PN method is combined with the discontinuous Galerkin (DG) finite elements method and yield to a full discretization of the multigroup neutron transport equation. The employed method is able to handle all geometries describing the fuel elements without any simplification nor homogenization. Moreover, the use of the matrix assembly-free method avoids building large sparse matrices, which enables producing high-order solutions in a small computational time and less storage usage. The resulting transport solver, called NYMO, has a wide range of applications; it can be used for a core calculation as well as for a precise 281-group lattice calculation accounting for anisotropic scattering. To assess the accuracy of this numerical scheme, it is applied to a three-dimensional (3-D) reactor core and fuel assembly calculations. These calculations point out that the proposed PN -DG method is capable of producing precise solutions, while the developed solver is able to handle complex 3-D core and assembly geometries.