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 discretizing 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 which can be done in parallel. In numerical examples we illustrate the algorithm's performance.