A novel numerical framework is discussed to simulate the time evolution of non-inertial particle size dis-tributions (or number density functions) in flames. The generic form of the population balance equa-tion is considered featuring nucleation, surface growth/loss and agglomeration/coagulation. This balance equation is first recast in a form that is prone to minimize spurious numerical errors in the simulation of surface growth/loss and collision integrals. Formally, this is achieved classifying the terms of the equa-tion into: (i) Lagrangian transport in size-space (surface growth/loss), (ii) relaxation rates of the parti-cle density at a given size (non-uniform growth/loss and negative contribution of collision integrals) and (iii) sources (nucleation and positive contribution of collision integrals). To secure accuracy, a high-order modal decomposition of the particle size distribution is introduced within every section of size consid-ered. A Legendre polynomials basis is used with Gauss-Lobatto quadrature points. By construction, the method performs very well for dealing with particle surface growth/loss and it is also highly accurate for the estimation of the collision integrals thanks to the high-order quadrature. This is confirmed simulat-ing canonical test cases of the literature to compare the numerical results against exact and analytical so-lutions. With a discretisation based on about 40 sections of size and with Legendre interpolation at the 5th-order, very good accuracy is obtained up to the third moment of the distributions for particle size ranging over up to 8 orders of magnitude. The method is cast to minimize computing cost. Strategies to couple this novel numerical framework with the simulation of carbon particles dynamics in flames are discussed.& COPY; 2022 The Combustion Institute. Published by Elsevier Inc. All rights reserved.