The probability distribution of neutron numbers in a symmetric subcritical reflected fissile sphere is numerically obtained using a one-speed diffusion approximation to the underlying backward Master equation. Employing an accurate space–time discretisation scheme, the coupled but closed system of equations for the number probabilities is sequentially solved as a function of position and time of an injected neutron. This solution is then used to construct the corresponding distributions for a random intrinsic source of arbitrary multiplicity. Numerical results clearly demonstrate the importance of including spatial dependence in the neutron number probability distributions, which show complex spatial behaviours that cannot be encapsulated in a point model system. This is especially evident in the case of the multi-region model where the presence of a reflector is seen to alter the approach to steady state of the number probabilities, while the material interface has a significant effect on the magnitude of the probabilities, both locally and globally.