A refinement of the analytic function expansion nodal (AFEN) method is described. By increasing the number of flux expansion terms in the way that the original basis functions are combined with the transverse-direction linear functions, the refined AFEN method can describe the flux shape in the nodes more accurately, since the added flux expansion terms still satisfy the diffusion equation. The additional nodal unknowns introduced are the interface flux moments, and the additional constraints required are provided by the continuity conditions of the interface flux moments and the interface current moments. Also presented is an algebraically exact method for removing the numerical singularity that can occur in any analytic nodal method when the core contains nearly no-net-leakage nodes. The refined AFEN method was tested on the Organization for Economic Cooperation and Development (OECD)-L336 mixed-oxide benchmark problem in rectangular geometry, and the VVER-440 benchmark problem and a nearly no-net-leakage node embedded core problem, both in hexagonal geometry. The results show that the method improves not only the accuracy in predicting the flux distribution but also the computing time, and that it can replace the corner-point fluxes with the interface flux moments without accuracy degradation, unless the problem consists of strongly dissimilar nodes. The possibility of excluding the corner-point fluxes increases the flexibility in implementing this method into the existing codes that do not have the corner-point flux scheme and may make it fit better for the nonlinear scheme based on two-node problems.