The following serial Fortran solvers for ordinary differential equation (ODE) initial value problems were written at LLNL. All are in the public domain and are freely available.
ODEPACK - A Systematized Collection of ODE Solvers
ODEPACK is a collection of Fortran solvers for the initial value problem for ordinary differential equation systems. It consists of nine solvers, namely a basic solver called LSODE and eight variants of it -- LSODES, LSODA, LSODAR, LSODPK, LSODKR, LSODI, LSOIBT, and LSODIS. The collection is suitable for both stiff and nonstiff systems. It includes solvers for systems given in explicit form, dy/dt = f(t,y), and also solvers for systems given in linearly implicit form, A(t,y) dy/dt = g(t,y). Two of the solvers use general sparse matrix solvers for the linear systems that arise. Two others use iterative (preconditioned Krylov) methods instead of direct methods for these linear systems. The most recent addition is LSODIS, which solves implicit problems with general sparse treatment of all matrices involved.
The ODEPACK solvers are written in standard Fortran 77, with a few exceptions, and with minimal machine dependencies. There are separate double and single precision versions of ODEPACK. The actual solver names are those given above with a prefix of D- or S- for the double or single precision version, respectively, i.e. DLSODE/SLSODE, etc. Each solver consists of a main driver subroutine having the same name as the solver and some number of subordinate routines. For each solver, there is also a demonstration program, which solves one or two simple problems in a somewhat self-checking manner.
Recently, the ODEPACK solvers were upgraded to improve their portability in numerous ways. Among the improvements are (a) renaming of routines and Common blocks to distinguish double and single precision versions, (b) use of generic intrinsic function names, (c) elimination of the Block Data subprogram, (d) use of a portable routine to set the unit roundoff, and (e) passing of quoted strings to the error message handler. In addition, the prologue and internal comments were reformatted, and use mixed upper/lower case. Numerous minor corrections and improvements were also made.
The above upgrade operations were applied to LSODE earlier than they were to the rest of ODEPACK, and the two upgrades were done somewhat independently. As a result, some differences will be apparent in the source files of LSODE and the other solvers -- primarily in the formatting of the comment line prologue of the main driver routine. In subroutines DLSODE/SLSODE and their subordinate routines, the prologue was written in "SLATEC style", while for the other solvers a more relaxed style was used. The differences are entirely cosmetic, however, and do not affect performance.
Documentation on the usage of each solver is provided in the initial block of comment lines in the source file, which (in most cases) includes a simple example. A demonstration program (in separate double/single precision versions) is also available.
A more detailed summary of ODEPACK, together with notes on usage documentation and installation, is supplied with the package. See also the papers ,  (for LSODE), and  (for LSODPK and LSODKR). (However, the document  does not reflect the upgrade operations described above.) For the convenience of users interested only in LSODE, that solver is being made available separately, as described below.
LSODE - SLSODE/DLSODE
LSODE (Livermore Solver for Ordinary Differential Equations) solves stiff and nonstiff systems of the form dy/dt = f(t,y). In the stiff case, it treats the Jacobian matrix df/dy as either a dense (full) or a banded matrix, and as either user-supplied or internally approximated by difference quotients. It uses Adams methods (predictor-corrector) in the nonstiff case, and Backward Differentiation Formula (BDF) methods (the Gear methods) in the stiff case. The linear systems that arise are solved by direct methods (LU factor/solve). LSODE supersedes the older GEAR and GEARB packages, and reflects a complete redesign of the user interface and internal organization, with some algorithmic improvements.
LSODE is available in separate double and single precision versions, called DLSODE and SLSODE. Documentation on the usage of DLSODE/SLSODE is provided in the initial block of comment lines in the source file, which includes a simple example. A demonstration program (in separate double/single precision versions) is also available.
VODE - SVODE/DVODE/ZVODE
VODE is a general purpose solver very similar to LSODE. However, it uses variable-coefficient methods (fixed-leading coefficient form) instead of the fixed-step-interpolate methods in LSODE. This and other features make it often more efficient than LSODE . VODE is available in separate double and single precision (real) versions, DVODE and SVODE, and also in a double precision complex version, ZVODE. For each solver -- DVODE, SVODE, ZVODE -- documentation on its usage is provided in the initial block of comment lines in the source file, which includes a simple example. A demonstration program is also available.
VODPK - DVODPK/SVODPK
VODPK is a variant of VODE which uses Preconditioned Krylov methods (which are iterative) for the linear systems that arise, instead of the direct methods in VODE. See  and . VODPK is available in separate double and single precision versions, DVODPK and SVODPK. Documentation on the usage of DVODPK/SVODPK is provided in the initial block of comment lines in the source file, which includes a simple example. A demonstration program is also available.
GEARBI is a variant of the older GEAR package (predecessor to LSODE). It solves stiff and nonstiff systems, using BDF and Adams methods. The unique feature of GEARBI is that, in the case of stiff systems, it uses a block-iterative method, Block-SOR, to solve the linear systems that arise at each time step. GEARBI is designed for use on problems that arise from the spatial discretization of PDE systems, such that the resulting ODE system has a regular block structure.
The original version of GEARBI, described and documented in , has been modified for the version here. The primary modification is in the treatment of the off-diagonal blocks of the Jacobian matrix. The new version assumes no structure for these blocks, but requires a user-supplied routine to compute matrix-vector products involving block-rows of the Jacobian, with the diagonal blocks excluded.
The GEARBI solver is provided in separate single and double precision versions. Both have a driver routine, called DRIVBI. Documentation on the usage of GEARBI is provided in the initial block of comment lines in the source file. An example program, with a sample output, is also supplied for each precision.
KRYSI is a solver for stiff systems, and is a variant of an implicit Runge-Kutta solver called SIMPLE. Both KRYSI and SIMPLE use the same 3-stage third order SDIRK method. But where SIMPLE uses a direct (dense) solver for the associated linear systems, KRYSI uses a preconditioned Krylov method (preconditioned GMRES iteration). See  for details.
The KRYSI solver is provided in separate single and double precision versions. Documentation on the usage of KRYSI is provided in the initial block of comment lines in the source file. An example program, with a sample output, is also supplied for each precision.