  We consider a linear inhomogeneous fractional evolution equation which is
obtained from a Cauchy problem by replacing its first-order time derivative
with Caputo's fractional derivative. The operator in the fractional evolution
equation is assumed to be sectorial. By using the inverse Laplace transform a
solution to the fractional evolution equation is obtained which can be written
as a convolution. Based on $L$-stable Runge-Kutta methods a convolution
quadrature is derived which allows a stable approximation of the solution.
Here, the convolution quadrature weights are represented as contour integrals.
On discretising these integrals, we are able to give an algorithm which
computes the solution after $N$ time steps with step size $h$ up to an
arbitrary accuracy $\varepsilon$. For this purpose the algorithm only requires
$\mathcal{O}(N)$ Runge-Kutta steps for a large number of scalar linear
inhomogeneous ordinary differential equations and the solutions of
$\mathcal{O}(\log(N)\log(\frac{1}{\varepsilon}))$ linear systems what can be
done in parallel. In numerical examples we illustrate the algorithm's
performance.
