A higher order nodal diffusion method is formulated, based on variational principle, Kantorovich’s variational method, and the patch test. In this framework, the relationship between finite element and nodal methods is discussed and the differences are pointed out. General, transverse integrated quasi-one-dimensional nodal equations are derived and matrix representation is given. In addition, a comparison with a similar approach is shown. A numerical solution is carried out using polynomial expansion of the source term and the corresponding analytic solution in alternating directions. Calculations of two-dimensional International Atomic Energy Agency and Biblis benchmark problems are performed and compared with results from the literature. It is shown that the first-order approximation yields the same order of accuracy as the standard nodal methods with quadratic leakage approximation, while the second-order approximation is considerably better.