Reaction networks in the bulk and on surfaces are widespread in physical, chemical and biological systems. In macroscopic systems, which include large populations of reactive species, stochastic fluctuations are negligible and the reaction rates can be evaluated using rate equations. However, with suitable couplings between them. In this paper we test the multiplane method and examine its applicability. We show that the method is accurate in the limit of small domains, where fluctuations are strong. It thus provides an efficient framework for the stochastic simulation of complex reaction networks with strong fluctuations, for which rate equations fail and direct integration of the master equation is infeasible. The method also applies in the case of large domains, where it converges to the rate equation results.