The population density approach has been used for modeling the dynamics of large-scale neuronal networks to consider the stochastic nature of the signal processing in the brain. From the conservation principle, the time evolution of population density can be described by a nonlinear partial differential equation (PDE). The finite difference method (FDM) has been proposed for solving this PDE. However, FDM is sensitive to the density gradient of the solution and is confined to problems that have a regular state space. An irregular state space is always obtained if realistic neuronal models are considered. In this study, the finite element method (FEM) is formulated to solve this PDE and apply it to solve the orientation tuning problem. The results show that when the state space discretization is coarse, FEM retains high accuracy whereas FDM does not. Furthermore, FEM reduces computation time by 90% compared to that required for FDM. In addition, FEM can easily handle the existence of a point source without any modifications. Due to its superior accuracy, efficiency, and consistency, FEM may be a better numerical technique for applying the population density approach to large-scale neuronal networks.