A numerical method for solving reaction rate equations

The numerical solution of differential equations for reacting species on a spatial grid is generally a difficult problem. The stiffness of the system of equations requires an implicit solver. Nonlinearities require that a quadratically convergent iterative method such as Newton-Raphson be used at each time step. This requires finding the Jacobian of the finite difference form of the differential equation, solving this matrix equation for each iteration, and storing this Jacobian in memory to allow vectorization of the problem. A method is described that eliminates these difficulties for such problems. The Newton iteration is replaced by the use of the Steffensen method, which is also quadratically convergent, but does not require the calculation of a Jacobian nor the solution of a matrix equation. This results in a simpler and faster solution algorithm, as well as less memory usage. The method is demonstrated by solving the rate equations for impurity species and electron temperature with ionization, recombination, ohmic heating and radiative power loss as source/sink terms. A description of a 2D resistive MHD code to employ this method is described