For several decades, microbubbles have been the primary choice for ultrasound contrast agents both for efficiency and safety reasons. To optimize their performance in various applications e.g. proton therapy, it is important to understand the dynamics of populations of nonlinearly oscillating microbubbles. Especially for high concentration clouds, multiple scattering can significantly influence the propagation of medical ultrasound. To numerically study the higher-order interaction between nonlinear microbubbles, we present a modified version of the Iterative Nonlinear Contrast Source (INCS) method. Based on a Neumann iterative scheme, the acoustic pressure is obtained by treating the nonlinearly scattering bubbles as individual contrast sources in a “linearized” background medium. In each iteration, the full nonlinear pressure field is updated by applying the 4D spatiotemporal convolution between the background Green’s function and the contrast sources obtained from the previous iteration. In this study, all microbubbles may exhibit an individual behavior. To accommodate this, in each iteration the solution of the Marmottant equation for each microbubble is obtained and used as the temporal signature of the respective contrast source. Using this scheme, each iteration adds an order of multiple scattering. The accuracy and the efficiency of INCS results will be demonstrated for monodisperse and polydisperse concentrations.