We consider simulation of multiple scattering of waves in sotropic and anisotropic media. The focus is on the construction of the phase function interpolation for the single scattering. The procedure is based on the construction of the adaptive partitioning of the angular variables that determine the phase function. The developed interpolation method allows us rather quickly to perform calculations for systems with very complicated phase function. Application of the proposed method is illustrated by calculating the multiple scattering of light in a nematic liquid crystal (NLC) which presents the uniaxial anisotropic system. For this system the grid corresponding to the adaptive partitioning is constructed and the transition to the diffusion regime for the photon distribution is presented.