<p><b>Abstract</b>—A nonlinear mathematical model is formulated for the two-step nucleation and growth of crystals in a supersaturated solution with allowance for the mass exchange between a crystallizing liquid and the environment. The mass exchange processes include the following two mechanisms: the removal of product crystals of a specified size from a supersaturated solution and the influx of impurities into the liquid to maintain its supersaturation. The integro-differential model of bulk crystallization consists of a first-order kinetic equation and an integral mass balance equation for a crystal radius distribution function and the solution supersaturation. The mathematical model also includes boundary and initial conditions. The growth rate of spherical particles, which is incorporated into the kinetic and balance equations, takes into account the two-step growth mechanism. Exact analytical solutions of the model are derived for a steady-state crystallization mode. Specifically, the steady-state solution supersaturation and the corresponding stationary crystal radius distribution function are determined. An approximate analytical solution for the nonstationary model is obtained using the Laplace transform in time and the saddle-point method. The solution is constructed in a parametric form, where the parameter is the maximum size of crystals grown in the solution by a certain time. The solution supersaturation dynamics is analyzed for various intensities of product crystal removal and impurity influx from outside. An increase in the crystal removal intensity is shown to cause a faster decrease in the liquid supersaturation, which is then restored due to the influx of impurities from outside. A nonstationary crystal radius distribution function and its moments are determined. The developed theory agrees with the experimental data.</p>