Abstract: Starting from the neutron transport equation, the diffusion equation is firstly derived via the P₁ approximation, and subsequently the neutron diffusion equation is obtained through the diffusion approximation. Compared with the diffusion equation, the P₁ equation directly retains the complete first-order anisotropic scattering and eliminates the diffusion approximation. Existing research indicates that for cases involving strong anisotropic scattering, the computational accuracy achieved using the P₁ equation is significantly higher than that obtained with the diffusion equation, without significantly increasing computation time. Therefore, to accurately solve this equation within the coarse-mesh nodal method framework, this paper referred to the variational nodal method for anisotropic scattering P_N equation, to derive the second-order even-parity form of the steady-state multigroup P₁ equation. Then, a theoretical model of the variational nodal method for solving the P₁ equation was constructed, which included functional formulation based on the Galerkin variational principle, Ritz discretization using three-dimensional orthogonal basis functions for the flux within nodes and for the current on surfaces, and variational derivation to yield the spatially discretized P₁ equation. Furthermore, based on an existing variational nodal method code for the diffusion equation, a variational nodal code for the P₁ equation, named NECP-VioletP1, was developed. The modifications included enhancing the cross-section input module to handle first-order scattering matrices, implementing new routines to calculate specific response matrices, and integrating the first-order scattering source term into the iterative solution scheme. And its’ verification and analysis were conducted through quantitative comparison with the finite difference code for solving the P₁ equation, set up for several benchmark cases including 2D single-assembly with different fuel enrichments and poison rod configurations, 2×2 assemblies with increasing heterogeneity, and 2D mini-core problem with a light water reflector. The results show that the maximum deviation of the effective multiplication factor (k_eff) calculated by this code from that of the finite difference code with refined-mesh is −12 pcm, and the maximum relative deviation of the fission rate distribution is −0.76%, with the largest errors predictably occurring at high-gradient regions like fuel-reflector interfaces. A comparative analysis on the mini-core problem highlights the significant advantage of the P₁ method over the diffusion method. Furthermore, spatial convergence studies indicate that a 5th-order flux expansion and 2nd-order current expansion within the variational nodal method achieve sufficient accuracy with computational efficiency superior to mesh refinement in the finite difference method (FDM). This study successfully establishes a robust and accurate variational nodal method for solving the P₁ equation. The developed NECP-VioletP1 code proves to be a highly effective solver, offering a superior accuracy compared to traditional diffusion theory. It provides a reliable core solver for advanced, high-fidelity core physics analysis based on the P₁ theory.