The process of particle nucleation and growth at the initial and intermediate stages of bulk crystallization in metastable liquids is studied. An integrodifferential model of balance and kinetic equations with corresponding boundary and initial conditions is formulated with allowance for non-stationary temperature/concentration field around each evolving particle. The model is solved using the saddle-point technique in a parametric form. The particle-radius distribution function, supercooling/supersaturation of liquid, total number of particles in liquid and their average size are found analytically. The melt supercolling (solution supersaturation) decreases with time due to the latent heat of phase transformation released by evolving crystals. As this takes place, the particle-radius distribution function is bounded by the maximal size of crystals and shifts to larger crystal radii with time as a result of particle nucleation and growth. © 2023, The Author(s), under exclusive licence to EDP Sciences, Springer-Verlag GmbH Germany, part of Springer Nature.