A new nodal method that directly solves the multidimensional diffusion equation without the transverse integration procedure is described. The new method expands the homogeneous flux distributions within a node in nonseparable analytic basis functions satisfying the neutron diffusion equation at any point of the node. Thus, the method accurately models large localized flux gradients in the vicinity of nodal corner points as well as nodal interfaces. To demonstrate its accuracy and applicability to realistic problems, the new method was tested on several well-known benchmark problems, including a mixed-oxide fuel problem, and the initial core of Ulchin Unit 1, which is a Framatome-type pressurized water reactor rated at 2775 MW (thermal). The results show that the new method significantly improves the accuracy in the nodal flux distribution and the core multiplication factor. The method also facilitates pin wise flux reconstruction since the homogeneous flux distributions obtained from the nodal calculation are very accurate and may be used directly in the reconstruction.