In this work, a cross-sectional averaged two-fluid model combined with a population balance model is applied to simulate the flow field and the bubble size distributions in a two-phase bubble column. The Martinez-Bazan breakage kernel and a modified Prince and Blanch coalescence kernel have been chosen to describe bubble breakage and bubble coalescence, respectively. In the present study, we discuss the Use of a higher-order spectral element method-the least-squares method-to compute the system of equations in a coupled manner. The least-squares method is highly accurate and has a number of advantages over the conventional numerical methods like the finite difference and finite volume methods. In contrast to the finite volume method, when desinging an overall solution algorithm, this least-squares method ensures that all the continuity equations are satisfied individually and it deals with both the convective and diffusive terms stably and accurately. The novel iterative algorithm solves the flow and the population balance equation ill a coupled manner. The model has been validated against experimental data obtained for two-phase flow in a bubble column. The predicted bubble size distribution and other flow quantities are in good agreement with the experimental data.