Mathematical models describing particle-size distributions in crystallization processes are of a challenging complexity. Depending on the considered physical phase system and the considered process, mathematical models of crystallization processes include phenomena such as primary nucleation, crystal growth, fines dissolution as well as agglomeration and/or breakage (attrition) of crystals. From a numerical point of view, agglomeration and breakage of crystals add particular (quadratic) complexity, since all particle sizes are connected by integral terms. A new powerful algorithm for the treatment of all these structures is presented in this paper. The method - called Galerkin h-p method - is based on a generalized finite-element scheme with self-adaptive grid- and order construction and is connected to a time discretization of Rothe's type. The algorithm can be applied to all combinations of the phenomena discussed above and needs no additional information on the form of the particle-size distribution. The numerical algorithm and the user-interface was implemented using object-oriented concepts leading to the simulation package PARSIVAL (PARTICLE SIZE EVALUATION). The paper presents some basic Features of PARSIVAL and the numerical algorithm as well as one illustrating example.