The prediction of the dispersed phase distribution plays a major role in multiphase chemical reactor engineering. The population balance equation (PBE) is a well-established equation for describing the evolution of the dispersed phase. However, the numerical solution of the PBE is computation intensive and challenging. In recent literature, the high-order least squares method has been applied to solve population balance (PB) problems. The interests in the least squares technique are based on the favorable numerical properties of the method. Moreover, by adopting a spectral method for the solution of the fundamental PBE, the statistical density function is directly obtained and the problem of reconstruction of the statistical density function is avoided, as is necessary, using moment methods. Furthermore, the least squares method is based on advanced linear algebra theory and thus is associated with involved implementation issues. For this reason, in this study, the theory and detailed implementation of the least squares method to solve the PBE for bubbly flow are outlined using an illustrated example.