In this article higher order in time numerical schemes with efficient time stepping for the solution of concentration and temperature distributions in fluidized beds using parallel computers are presented. The mathematical model equations consist of strongly coupled and semi linear convection-diffusion-reaction equations. Invariant regions for the model are derived to check the solution bounds. The numerical discretization for the space using the finite element method is presented and the numerical treatment is enhanced by using adaptive and higher order linearly implicit Runge-Kutta methods for the time discretization. For different time stepping methods and different spatial grid sizes numerical results are obtained and compared. The methods used show a clear improvement for the problem under consideration compared to previously presented results (Nagaiah, Warnecke, Heinrich, & Peglow, 2008). Additionally, the higher order time stepping methods yield a good parallel efficiency, paving the way for the efficient study of more complex phenomena. (C) 2012 Elsevier Ltd. All rights reserved.