The time and phase-space dependent backward master equation is used to develop and numerically solve a coupled system of transport equations for the probability distribution of the neutron number in subregions of a spherically symmetric, reflected, subcritical plutonium sphere. The number distributions are computed for a single initial neutron injected into the assembly and localized in phase space as well as in the presence of a uniformly distributed spontaneous fission source in the fissile region. A standard multigroup, discrete ordinates scheme with second-order spatial and fully implicit time discretization proved sufficiently accurate for this application. The results presented show complex behaviors arising from the material interface and spectral effects due to neutron slowing down that cannot be encapsulated in a lumped model. Additionally, low-order spatial moments were computed both by averaging the number distributions of finite order and directly solving the transport equations for the moments using the same numerical scheme. While generally excellent agreement is observed between the two approaches, the truncation order has a noticeable effect on the accuracy of the higher moments that are computed using the number distributions.