  Over the billions of years since the Big Bang, the lives, deaths and
afterlives of stars have enriched the Universe in the heavy elements that make
up so much of ourselves and our world. This review summarizes the methods used
to evolve these nuclear abundances within astrophysical simulations. These
methods fall into 2 categories; evolution via rate equations and via
equilibria. Because the rate equations in nucleosynthetic applications involve
a wide range of timescales, implicit methods have proven mandatory, leading to
the need to solve matrix equations. Efforts to improve the performance of such
rate equation methods are focused on efficient solution of these matrix
equations, in particular by making best use of the sparseness of these
matrices, and finding methods that require less frequent matrix solutions.
Recent work to produce hybrid schemes which use local equilibria to reduce the
computational cost of the rate equations is also discussed. Such schemes offer
significant improvements in the speed of reaction networks and are accurate
under circumstances where calculations which assume complete equilibrium fail.
