Reliable and accurate simulations of the thermal runaway of Li-ion batteries are crucial for designing thermal management system. However, the stiff nature of thermal abuse kinetics makes it computationally challenging to obtain numerical solutions. To address this, we employed modified Patankar–Runge–Kutta (MPRK) methods for thermal runaway simulations. The explicit nature of these methods leads to computationally efficient simulations, while their ability to maintain positive physical quantities ensures robust integrations. We first investigated the performance of the second- and third-order MPRK methods, along with the well-known fourth-order Runge–Kutta method, for zero-dimensional thermal runaway simulations. Numerical tests revealed the superiority of the third-order MPRK method in terms of robustness and accuracy. To further improve its accuracy during a rapid temperature rise, we introduced an adaptive time-stepping method. This adaptive third-order MPRK method was then implemented in Ansys Fluent for three-dimensional simulations of thermal runaway propagation. Benchmark tests demonstrated an increase of 10 times in the computational speed compared with Fluent’s built-in solver, which was achieved by allowing significantly larger thermal time steps without compromising accuracy. The robust and efficient computational method proposed in this study will pave the way for future studies on multiphysics simulations of thermal runaway phenomena.