Coupled sets of ordinary differential equations (ODEs) are used to describe the evolution of the interactions among chemical species as well as many other local processes. ODEs appear, for example, when spectral and other expansion methods are used to solve timedependent partial differential equations. In these cases, spatial derivatives are converted to algebraic relationships leaving ODEs to be integrated in time. ODEs also describe the motions of projectiles and orbiting bodies, population dynamics, electrical circuits, local temperature equilibration, momentum interchange among phases in multiphase flows, the decomposition of radioactive material, and energy level and species conversion processes in atomic, molecular, and nuclear physics.

Algorithms for integrating ODEs were not originally derived by numerical analysts or applied mathematicians, but by scientists interested in solving specific sets of equations for their particular applications. Bashforth and Adams (1883), for example, developed a method for solving the equations describing capillary action. One of the first algorithms to cope with the difficulties of integrating stiff ODEs was suggested by Curtiss and Hirschfelder (1952) for chemical kinetics studies. Ten years after Curtiss and Hirschfelder identified the stiffness problem in ODEs, Dahlquist (1963) identified numerical instability as the cause of the difficulty and provided basic definitions and concepts that are still helpful in classifying and evaluating algorithms. The importance of the practical applications has spurred active research in developing and testing integration methods for solving coupled ODEs. Continued efforts of applied mathematicians have put the numerical solution of ODEs on a sounder theoretical basis and have provided insights into the constraints imposed by stability, convergence, and accuracy requirements.