  The finite difference scheme with the shifted Gr\"{u}nwarld formula is
employed to semi-discrete the fractional diffusion equations. This spatial
discretization can reduce to the large system of ordinary differential
equations (ODEs) with initial values. Recently, boundary value method (BVM) was
developed as a popular algorithm for solving large systems of ODEs. This method
requires the solutions of one or more nonsymmetric, large and sparse linear
systems. In this paper, the GMRES method with the block circulant
preconditioner is proposed for solving these linear systems. One of the main
results is that if an $A_{\nu_1,\nu_2}$-stable boundary value method is used
for an m-by-m system of ODEs, then the preconditioner is invertible and the
preconditioned matrix can be decomposed as I+L, where I is the identity matrix
and the rank of L is at most $2m(\nu_1+\nu_2)$. It means that when the GMRES
method is applied to solve the preconditioned linear systems, the method will
converge in at most $2m(\nu_1+\nu_2)+1$ iterations.Finally, extensive numerical
experiments are reported to illustrate the effectiveness of our methods for
solving the fractional diffusion equations.
