Green’s function nodal expansion method (GNEM) is developed for the efficient numerical solution of the multidimensional neutron diffusion equation. It is an improved version of the nodal expansion method (NEM) and the nodal Green’s function method (NGFM). The node interior fluxes are approximated by a high-order polynomial expansion as in NEM. The nodal surface fluxes are coupled with the net currents by using the Green’s function method to improve accuracy. The GNEM computer code is encoded and tested. The numerical results demonstrate that GNEM has the same accuracy as NGFM while it is twice as fast as NGFM.