In this paper, we are concerned with numerical methods for the solution of initial-boundary value problems of anomalous diffusion equations of order α. ∈. (1, 2). The classical Crank-Nicholson method is used to discretize the fractional diffusion equation and then the spatial extrapolation is used to obtain temporally and spatially second-order accurate numerical estimates. Two preconditioned iterative methods, namely, the preconditioned generalized minimal residual (preconditioned GMRES) method and the preconditioned conjugate gradient for normal residual (preconditioned CGNR) method, are proposed to solve relevant linear systems. Numerical experiments are given to illustrate the efficiency of the methods.