A robust, efficient numerical method for computing the effectiveness factor of a heterogeneous reaction in a catalyst is developed in this study. The method is based on shooting at the outer surface of the catalyst and is optimized for an accurate estimation of the concentration gradient at the outer surface. The shooting at the outer surface, however, is inherently unstable for a cylindrical or a spherical catalyst, and it also can become unstable with an increasing Thiele modulus even for a slab catalyst. From an analysis of the governing equation, however, three criteria are developed to make the method work efficiently in the presence of instability. The numerical method is shown to be effective for diverse kinetic expressions, from simple power-law to sophisticated Langmiur-Hinselwood kinetics, and for isothermal or non-isothermal catalysts. The method is also shown to be easily applicable to a more complex case of multiple steady states, estimating all the corresponding effectiveness factors. In this method, a simple rule that determines the stability of each steady state is proposed. (C) 2004 Elsevier Ltd. All rights reserved.