A nodal method for the solution of the multidimensional neutron diffusion equation is developed and evaluated. The method is based on the linear form of the nodal balance equation written in terms of the average partial currents across the surfaces of the node. Green's functions for one-dimensional in-group diffusion-removal operators are used to generate a coupled set of one-dimensional integral equations defined over a subdomain or node. These integral equations represent an exact (local) solution to the coupled set of one-dimensional differential equations obtained by spatially integrating the multidimensional diffusion equation over directions transverse to each coordinate direction. The integral equations are approximated using a weighted residual procedure applied within each node. The resulting matrix equations, when solved in conjunction with the linear form of the nodal balance equation, provide the necessary additional relationships between the interface partial currents and the flux within the node. The nodal method is applied to several two- and three-dimensional light water reactor benchmark problems and to a four-group liquid-metal fast breeder reactor problem. These results demonstrate the capability of the method to yield very accurate steady-state and transient results in significantly smaller computing times than those required by standard finite difference methods.