The authors study the discretization in time of initial-value problems for \[ {\partial u\over\partial t}+ \int^t_0 \beta(t- s)Au(s)\,ds= f(t),\quad t> 0, \] where \(\beta(t)= t^{\alpha-1}/\Gamma(\alpha)\) \((0< \alpha< 1)\) and \(A\) represents an elliptic differential operator in space. Based on the approach introduced by D. Sheen, I. H. Sloan and V. Thomée [Math. Comput. 69, No. 229, 177–195 (2000; Zbl 0936.65109) and IMA J. Numer. Anal. 23, No. 2, 269–299 (2003; Zbl 1022.65108)], the solution is represented in the form of an integral along a smooth curve extending into the half-plane \(\text{Re}(z)< 0\) of \(\mathbb{C}\) and the integral then is approximated by appropriate quadrature formulas (given by two different adaptations of the trapezoidal, rule). The resulting systems of complex elliptic equations can be solved in parallel. Following the derivation of high-order error estimates (for the two different choices of quadrature rules) and a discussion of finite-element spatial discretization, the theoretical results are illustrated by two numerical examples.