An adaptive orthogonal collocation on finite elements method with adaptively varied upper bound of the relevant size interval is developed for numerical solution of population balance equation of crystallizers. Nucleation producing monosized and heterosized nuclei, size-dependent crystal growth, seeding, classified product removal and fines removal with dissolution are included into the model. Adaptation of the number and length, as well as the distribution of finite elements over the variable length computational interval is carried out forming a number of adaptation rules, based on ordering the finite elements of size coordinate according to the maxima of first derivatives of the population density function. The approximation is obtained using the Lagrange interpolation polynomials. The method is used for solving the mixed set of nonlinear ordinary and partial differential equations, forming a detailed dynamical model of continuous crystallizers with product classification and/or fines removal. The program can be used efficiently for simulation of stationary and dynamic processes of crystallization systems, computing either transients or long-time oscillating steady states generated by different nonlinear phenomena, and internal and external feedbacks. (c) 2007 American Institute of Chemical Engineers AIChE J, 53: 3089-3107, 2007.