The popular backward differentiation methods for solving stiff systems of ordinary differential equations, being implicit, require the solution of a linear algebraic system at each time step. The coefficient matrix is closely related to the Jacobian matrix of the differential system, and large systems may require considerable storage for the Jacobian. The authors consider an alternative to the usual direct methods for solving the algebraic system which does not require the storage of the coefficient matrix in any form. The method is an iterative one, in which a Krylov-subspace projection method known as the incomplete orthogonalization method is combined with the usual Newton iteration, giving an approach in which the linear system in the Newton iteration is solved only approximately. The authors analyse this class of methods and present algorithms for its incorporation into the LSODE package. The numerical examples which they report indicate that the new approach can be advantageous, though further research is required for the class of problems for which the method is best suited.