A novel numerical framework is discussed to simulate the time evolution of non-inertial particle size distributions (or number density functions) in flames. The generic form of the population balance equation 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 equation into: (i) Lagrangian transport in size-space (surface growth/loss), (ii) relaxation rates of the particle 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 considered. 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 simulating canonical test cases of the literature to compare the numerical results against exact and analytical solutions. 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.