Summary: We present an extrapolation type algorithm for the numerical solution of fractional order differential equations. It is based on the new result that the sequence of approximate solutions of these equations, computed by means of a recently published algorithm by K. Diethelm
[ETNA, Electron. Trans. Numer. Anal. 5, 1-6 (1997; Zbl 0890.65071
)], possesses an asymptotic expansion with respect to the stepsize. From this we conclude that the application of extrapolation is justified, and we obtain a very efficient differential equation solver with practically no additional numerical costs. This is also illustrated by a number of numerical examples.