A three-dimensional variational nodal method (VNM) is presented for pressurized water reactor core calculations without fuel-moderator homogenization. The nodal functional is presented and discretized to obtain response matrix equations. Within the nodes, finite elements in the x-y plane and orthogonal polynomials in z are used to approximate the spatial flux distribution. On the lateral interfaces, orthogonal polynomials are employed. On the axial interfaces, the finite elements facilitate a spatially accurate current representation that has proven to be a challenge for the method of characteristics–based two-dimensional/one-dimensional approximations which typically rely on spatial homogenization. The angular discretization utilizes an even-parity integral method within the nodes, with the integrals evaluated using high-order Chebyshev-Legendre cubature. On the lateral and axial interfaces, low-order spherical harmonics (Pn) are augmented by high-order Pn expansions to which quasi-reflected conditions are applied. With quasi-reflected conditions, the solution converges to the high-order Pn solution for an infinite lattice of identical cells with no gradient, while the low-order Pn expansions handle global gradients in both the radial and axial directions. The method is implemented in the PANX code and applied first to a number of model problems to study convergence of the space-angle approximations and then to the C5G7 benchmark problems. Multigroup Monte Carlo solutions provide reference values for eigenvalues and pin-power distributions.