Here we present results of Monte Carlo simulations employing an effective numerical approach, based upon a description of radiative transfer in terms of the Bethe-Salpeter equation. We develop a simulation algorithm which presents a stochastic calculation of the multi-fold spatial integrals of the Bethe-Salpeter iterative solution. The algorithm is a modification of the widely known MCLM method with a signal registration recipe resulting in a significant reducing of sampling volume. For isotropic scattering the obtained results exhibit an excellent agreement with the known theoretical solution. We simulate backscattering from layers of different thickness for the Henyey-Greenstein or Rayleigh-Gans scattering patterns.