The nonlocal nature of the fractional integral makes the numerical treatment of fractional dierential equations expensive in terms of computational eort and memory requirements. In this paper we propose a method to reduce these costs while controlling the accuracy of the scheme. This is achieved by splitting the fractional integral of a function f into a local term and a history term. Observing that the history term is a convolution of the history of f and a regular kernel, we derive a multipole approximation to the Laplace transform of the kernel. This enables the history term to be replaced by a linear combination of auxiliary variables dened as solutions to standard ordinary dierential equations. We derive a priori error estimates, uniform in f, and obtain estimates on the number of auxiliary variables required to satisfy an error tolerance. The resulting formulation is discretized to produce a time stepping method. The method is applied to some test cases to illustrate the performance of the scheme.