A highly efficient numerical scheme has been incorporated to solve a traditional two-fluid model for dense gas-solid flow in large scale fluidized beds. Difficulties associated with the numerical solution and boundary condition enforcement especially in the azimuthally direction and at the axis of the cylindrical domain are discussed in detail. Higher order discretization schemes with deferred correction approach have been implemented for the convection terms. The p-epsilon algorithm (Van der Hoef er al., 2006) has been implemented to deal with densely packed regions in the fluidized bed. A modified SIMPLE algorithm is used to solve the pressure and volume fraction corrections, and a projection method is used to obtain solutions of the momentum equations. The resulting semi-implicit method allows for using larger time steps and produces accurate results in a stable and efficient manner. Numerical tests on bubbling fluidized beds are undertaken and compared with experimental data reported by Laverman et al, 2012. The simulation results are found to be in good agreement. In addition a comparative study has been performed quantifying the effects of grid size, flux limiters, frictional model and coefficient of restitution. (C) 2013 Elsevier Ltd. All rights reserved.