The R-function method is applied to model the geometry of the critical facility ZED-2 as an illustrative example. Each material region is represented by a real, continuous and differentiable function of spatial coordinates that is referred to as the domain function. Instead of using particular functions for individual elements, a single function is constructed to represent all spatial domains that contain the same material, for instance all fuel pins, etc. Owing to continuity and differentiability, the domain functions can be used to construct a basis for approximate solution of the related boundary-value problems. A study of the approximation ability of the domain functions is carried out using the two-group spatial neutron flux distribution in the ZED-2 facility as a model problem. The least squares method is used to determine the unknown coefficients by minimizing the discrepancy between the reference MCNP solution and a power series of domain functions. The results show that a modest number of terms in the series is able to produce a good approximation of the neutron flux distribution specified on a mesh grid of 500 x 500 points.