Various solution methods were studied for the solution of a general isothermal dynamic bubble column reactor model. The general dynamic reactor model can be used to describe various degrees of backmixing in co- and countercurrent bubble columns. The solution methods were all based on discretization of the transient film-penetration model in the film dimension and on discretization of the column axial dimension. The resulting ODE-system was then integrated with a stiff integration algorithm applying sparse matrix techniques. The solution methods were tested using a single first order irreversible reaction and its autocatalytic modification. Furthermore, more complicated reaction kinetics of two competitive-consecutive second order reactions referring to the chlorination of para-cresol was simulated. The results show that methods based on orthogonal collocation perform extremely well in the film dimension, but suffer from serious stability problems in the column dimension. The discretization methods based on a finite difference approach are robust and reliable, but numerous grid points must be applied in order to obtain an acceptable accuracy. The multipoint finite difference formulae can be used to enhance the accuracy considerably. Also, a dynamic formulation of a semi-empirical outlet boundary condition for the axial dispersion model was tested. The proposed formulation is numerically superior to the traditional Danckwerts’ outlet boundary condition in dynamic column simulation with all tested kinetics.