We study the solution of neutral delay differential equations (ND-DEs) by using boundary value methods (BVMs). The BVMs require the solution of nonsymmetric, large and sparse linear systems. The GMRES method with Strang-type block-circulant preconditioner is proposed to solve these linear systems. We show that, if an A-stable BVM is used for solving a system of NDDEs, then our preconditioner is invertible and the spectrum of the preconditioned system is clustered. It follows that, when the GMRES method is applied to the preconditioned systems, the method can converge rapidly. Numerical results are given to show that our method is effective.