A new numerical approach for solving population balance equations (PBE) is proposed and validated. The method employs a combination of basis functions, defined on finite elements, to approximate the sought distribution function. Similarly to other methods of the same family, the PBE are solved only in a finite number of values of the internal coordinate (grid points). The peculiarity of the method is the use of a logarithmic, shape-preserving interpolation (LSPI) procedure to estimate the values of the distribution in between grid points. The main advantages of the LSPI method compared to other approaches of the same category are: (i) the stability of the numerical approach (i.e., the absence of oscillations in the distribution function occurring when using "standard" cubic splines and a low number of elements), and (ii) the conceptual and implementation simplicity, as no mathematical manipulation of the PBE is required. (C) 2018 Elsevier Ltd. All rights reserved.