  Aims: The aim of this work is to compare the performance of three reaction
network integration methods used in stellar nucleosynthesis calculations. These
are the Gear's backward differentiation method, Wagoner's method (a 2nd-order
Runge-Kutta method), and the Bader-Deuflehard semi-implicit multi-step method.
  Methods: To investigate the efficiency of each of the integration methods
considered here, a test suite of temperature and density versus time profiles
is used. This suite provides a range of situations ranging from constant
temperature and density to the dramatically varying conditions present in white
dwarf mergers, novae, and x-ray bursts. Some of these profiles are obtained
separately from full hydrodynamic calculations. The integration efficiencies
are investigated with respect to input parameters that constrain the desired
accuracy and precision.
  Results: Gear's backward differentiation method is found to improve accuracy,
performance, and stability in integrating nuclear reaction networks. For
temperature-density profiles that vary strongly with time, it is found to
outperform the Bader-Deuflehard method (although that method is very powerful
for more smoothly varying profiles). Wagoner's method, while relatively fast
for many scenarios, exhibits hard-to-predict inaccuracies for some choices of
integration parameters owing to its lack of error estimations.
