Population balance models describing particle formation kinetics in supercritical crystallization processes in general and of gas antisolvent crystallization (GAS) process in particular are rather complicated. Depending on the considered particle formation mechanism, population balance model of GAS process may include phenomena such as primary nucleation, secondary nucleation, crystal growth, as well as agglomeration and/or breakage (attrition) of crystals. The resulting equation is often of the integro-partial differential form. A powerful numerical algorithm for the treatment of the implemented population balance model structures is presented in this paper. The method/algorithm is a combination of the Lax-Wendroff and Crank-Nicholson methods. Algorithm simulations were performed for changes in the main GAS process operating parameters, i.e., the antisolvent addition rate and saturation level. The simulations were performed at a process temperature of 25 degrees C, while the antisolvent addition rate, Q(A), was varied between 1 and 100 ml/min, and the initial solute concentration was varied between 25% and 100% of the concentration ratio. The numerical algorithm was successfully able to treat the implemented population balance model structures, producing smooth dynamic and steady state particle size distributions. Moreover, simulation findings were reasonably consistent with the experimentally obtained data. (C) 2009 Elsevier Ltd. All rights reserved.