The `multiphase particle-in-cell coupled with population balance equation' (MP-PIC-PBE) method is introduced for simulating multi-scale multiphase particulate flows. This method couples the meso-scale fluid dynamics simulated by the MP-PIC method with the simulation of the micro-scale particle size distribution. The homogeneous population balance equation is calculated for each discrete particle tracked in a Lagrangian frame, after the MP-PIC numerical procedure is followed at each time instance. This approach allows the particulate phase to accommodate the particulate stresses using spatial gradients and allows the Lagrangian description to predict particle properties by the PBE. For the antisolvent crystallization of Lovastatin in a biradial mixer, the proposed method is compared to an existing method that simulates the spatiotemporal evolution of the particle distribution by combining a multi-environment probability density function with the spatially varying PBE. The MP-PIC-PBE method has lower computational cost and provides more detailed information, such as particle age and location.