A two-stage third-order numerical scheme is proposed for solving ordinary differential equations. The scheme is explicit and implicit type in two stages. First, the stability region of the scheme is found when it is applied to the linear equation. Further, the stability conditions of the scheme are found using a linearized homogenous set of differential equations. This set of equations is obtained by applying transformations on the governing equations of heat and mass transfer of incompressible, laminar, steady, two-dimensional, and non-Newtonian power-law fluid flows over a stretching sheet with effects of thermal radiations and chemical reaction. The proposed scheme with an iterative method is employed in two different forms called linearized and non-linearized. But it is found that the non-linearized approach performs better than the linearized one when residuals are compared through plots. Additionally, the proposed scheme is compared to the second-order central finite difference method for second-order non-linear differential equations and the Keller-Box/trapezoidal method for a linear differential equation. It is determined that the proposed scheme is more effective and computationally less expensive than the standard/classical finite difference methods. Moreover, the impact of magnetic parameter, radiation parameter, modified Prandtl and Schmidt numbers for power-law fluid, and chemical reaction rate parameter on velocity, temperature, and concentration profiles are displayed through graphs and discussed. The power-law fluid’s heat and mass transfer simulations are also carried out with varying flow behavior index, sheet velocity, and mass diffusivity. We hoped that this effort would serve as a guide for investigators tasked with resolving unresolved issues in the field of enclosures used in industry and engineering.