Evolution of crystal ensembles in supersaturated solutions is studied at the initial and intermediate stages of bulk crystallization. An integro-differential model includes fluctuations in crystal growth rates, initial crystal-size distribution and arbitrary nucleation and growth kinetics of crystals. Two methods based on variables separation and saddle-point technique for constructing a complete analytical solution to this model are considered. Exact parametric solutions based on these methods are derived. Desupersaturation dynamics is in good agreement with the experimental data for bovine and porcine insulin. The method based on variables separation has a strong physical limitation on exponentially decaying initial distribution and leads to the distribution function increasing with time. The method based on saddle-point technique leads to a dome-shaped crystal-size distribution function decreasing with time and has no strong physical limitations. The latter circumstance makes this method more reasonable for describing the kinetics of bulk crystallization in solutions and melts.