We consider the solution of a system of ordinary differential equations (ODEs) by waveform relaxation (WR) iterations in conjunction with boundary value methods (BVMs). The WR method is a continuous-in-time analogue of the stationary method and it iterates with functions. In each WR iteration, we use BVMs to discretize systems of ODEs. BVMs are relatively new ODE solvers based on linear multistep formulae. In this paper, we discuss the use of the generalized minimal residual (GMRES) method with block-circulant-circulant-block preconditioners for solving the linear systems arising from the application of BVMs in each WR iteration. These preconditioners are effective in speeding up the convergence rate of the GMRES method. Numerical experiments are presented to illustrate the effectiveness of our methods.