For solving a system of stiff ordinary differential equations numerically, implicit integrators with infinite stability regions are likely to work with relatively large step sizes when the accuracy requirements are not too stringent. A systematic way of exploiting sparsity when integrating atmospheric chemistry equations is developed. The authors discuss and evaluate reordering techniques that lead to minimal fill-in during LU decomposition. They then test various linear system solvers, in particular showing that the chosen routine is twice as fast as the one used by J. Verwer, J. G. Blom, M. van Loon, and E. J. Spee [Atmos. Environ. 30, 49 (1996)]. These ideas are tested on two comprehensive chemical mechanisms used to study stratospheric and tropospheric chemistry.