Abstract Multiphase flows play an important role in the chemical and process industries and significant strides have recently been achieved in the design of such systems using the population balance models. However, some uncertainties still remain concerning the stability and accuracy of the numerical solution of such integro-differential equations. This paper proposes a new methodology for solving the discretized population balance equation by minimizing the finite domain errors that often arise while discretizing the drop size domain. It relies on the use of the size distribution sampling approach combined with a moving grid technique. In addition, an enhanced solution stability algorithm was proposed and which relies on monitoring the onset of errors in the various birth and death terms encountered in PBE. It consequently allows for corrective action to be undertaken before the errors propagate in an uncontrollable fashion, and was found to improve the stability and robustness of the solution method even under very high shear rate conditions. The proposed algorithm was tested using the model of Coulaloglou and Tavlarides (1977) under breakage and coalescence dominated conditions in low, moderate and high energy dissipation regions, and was found to provide a stable solution that accurately predicts the quasi-equilibrium Sauter mean diameter.