In this thesis, we study the solution of linear systems arising from numerical methods for solving ordinary differential equations (ODEs). We use boundary value methods (BVMs) to discretize ODEs. These implicit numerical methods for solving ODEs require the solutions of non-symmetric, large and sparse linear systems. Therefore the preconditioned iterative methods were proposed instead of the direct methods. We consider the generalized minimal residual (GMRES) method, one kind of Krylov subspace methods, and the waveform relaxation (WR) method in this thesis. Some basic theories of these methods and circulant-like preconditioners are given in Chapter 1. In Chapter 2, we study the use of the GMRES method with circulant-like preconditioners to solve the linear systems generated from the discretization of ODEs. Particularly, a new type of circulant-like preconditioners, circulant- block (CB) preconditioners, is proposed for such problems. We show that the operation cost is lower than those of the block-circulant (BC) preconditioners and the block-circulant-circulant-block (BCCB) preconditioners in some cases. The estimate of the convergence rate is also given. In addition, numerical examples are presented to illustrate the effectiveness of this method. In Chapter 3, we consider another iterative method, the WR method, for solving ODEs. Since in each WR iteration, we are required again to solve a system of ODEs, we apply the GMRES method with Strang's BC preconditioner to solve the linear systems arising from the application of BVMs to each WR iteration. We briefly study the error estimate of WR iterations. Under stability assumption on a given BVM, we then prove that the Strang's preconditioner is invertible. In this chapter, we also study the operation cost and convergence rate of the preconditioned GMRES methods with respect to different types of WR versions. Numerical results are also given.