Abstract
We are concerned with an efficient numerical solution of linear equations at each time-stepping of the trapezoidal rule applied to a system of linear ordinary differential equations (ODEs) with a constant coefficient matrix of large dimension. We do not assume that the matrix is symmetric. Hence numerical solutions in the family of BiCG method are sought. The present paper describes a method to reuse Krylov subspaces in the BiCGSTAB process over a number of computational steps. It can suppress increase of the memory usage as well as reduce the total number of BiCGSTAB iterations. Numerical examples depict its efficiency.