Transport and reaction phenomena in catalytic pellet reactors are often difficult to analyze because of coupling between heat and mass transport occurring at different space and time scales. To calculate the reactor concentrations and temperatures, it is necessary to account for the species reaction and transport occurring in the reactor bulk at the macroscopic level as well as the catalyst pellets at the microscopic level. The resulting approach yields a large system of nonlinear partial differential equations with multiple scales and solutions that are difficult to find numerically. In addition, the catalyst pellets may operate in multiple steady states for identical conditions. Conventional computational methods may entirely miss the multiplicity phenomenon at the catalyst pellet level and, as a result, may not correctly predict overall reactor yields. In this paper, we introduce two numerical techniques to address multiple scales and multiplicity in heterogeneous reaction models. The first method expands existing bisection with "shooting"; the second global method deploys orthogonal collocation over finite elements with niche evolutionary algorithms. We also propose a new multiscale method entitled effectiveness factor maps to expedite and simplify the numerical effort to solve transport and reaction phenomena at different length scales.