The authors investigate both theoretical and practical advantages of incorporating codes for the Krylov subspace projection method, for solution of systems of linear equations, into implicit codes for large systems of stiff ordinary differential equations. At each mesh point an implicit code requires solution of a linear system of equations whose coefficient matrix typically has the form \(A=I-h\beta_ 0J\), \(\beta_ 0=const.\), where J is the Jacobian matrix associated with the differential system. At points where the differential system is stiff the dominant part of the solution of the linear system will lie in the eigenspaces of A associated with those eigenvalues having greatest real part.
The Krylov subspace method, which employs projections onto the eigenspaces of A, is well suited to solution of linear systems of this sort. At points where the differential system is not stiff alternatives to the Krylov algorithm may be employed. Assuming the eigenvalues of A to be real and distinct the authors present error bounds for the Krylov method and they give sufficient conditions for convergence of an iterative procedure based on this method. Results of trial computations are discussed in which the basic code used for the differential system was LSODE modified to employ a combination of the Krylov subspace methods of W. E. Arnoldi [Quart. Appl. Math. 9, 17-29 (1951; Zbl 0042.128)] and of the second author [SIAM J. Numer. Anal. 19, 485-506 (1982; Zbl 0483.65022)]. The results are reported to be favorable in terms of running time and core memory requirements when they are compared with computations using LSODE with the Gaussian elimination algorithm replacing the Krylov subspace algorithm.