The distribution of photons is described by the integrodifferential Boltzmann equation. However, this equation is difficult to solve because of the scatter integral and the size of the problem when discretized. A method is presented that models all orders simultaneously using spherical harmonics to efficiently describe the directional dependence. The method also partitions the problem to take advantage of the fact that a scattering event results in an energy loss. Then, within a parallel computing paradigm, it is assumed that the source and the attenuation map of the medium are known, and an attempt is made to solve for the distribution. This distribution can then be used in any imaging system model to approximate experimental images.