The R-function method is applied to the multidimensional steady-state neutron diffusion equation. Using a variational principle, the nested element approximation is formulated. Trial functions taking into account the geometric shape of material regions are constructed. The influence of both the surrounding regions and the corner singularities at the external boundary is incorporated into the approximate solution. Benchmark calculations show that such an approximation can yield satisfactory results. Moreover, in the case of complex geometry, the presented approach would result in a significant reduction in the number of unknowns compared to other methods.