A spatially and temporally adaptive completely implicit finite difference algorithm, which can accurately capture the moving temperature and reaction fronts in fixed-bed reactors with highly exothermic reactions, is developed and presented. For spatial adaptation linear interpolation proved to be more robust than higher-order interpolations. Node placement based on the magnitude of the second spatial derivative was implemented and shown to lead to more accurate results and a more robust algorithm than the traditional node placement based on the magnitude of the first derivative. For improved accuracy the interdependence of the spatial mesh size and the time step size is taken into account by coupling spatial and temporal adaptation. Choice of the time step is also related to the physics of the problem. Advantages of the developed algorithm over traditional finite difference and polynomial approximation methods are discussed. The validity of the developed algorithm is proven by comparing the computed results with analytical solution of a linear heat regenerator model. A nonlinear model of a fixed-bed reactor with highly exothermic reaction occurring in the gas phase is also solved and the results are presented. The reasonableness of the solution for this nonlinear model is shown by comparison of model calculated and approximate analytical equation for the front velocity. The robustness of the developed algorithm in solving multiple coupled equations is demonstrated by solving a nonlinear model of a fixed-bed reactor in which both exothermic and endothermic reactions occur in the gas phase.