Abstract An initial-boundary value problem, whose differential equation contains a sum of fractional time derivatives with orders between 0 and 1, is considered. Its spatial domain is ( 0 , 1 ) d {(0,1)^{d}} for some d ∈ { 1 , 2 , 3 } {d\in\{1,2,3\}} . This problem is a generalisation of the problem considered by Stynes, O’Riordan and Gracia in SIAM J. Numer. Anal. 55 (2017), pp. 1057–1079, where d = 1 {d=1} and only one fractional time derivative was present. A priori bounds on the derivatives of the unknown solution are derived. A finite difference method, using the well-known L1 scheme for the discretisation of each temporal fractional derivative and classical finite differences for the spatial discretisation, is constructed on a mesh that is uniform in space and arbitrarily graded in time. Stability and consistency of the method and a sharp convergence result are proved; hence it is clear how to choose the temporal mesh grading in a optimal way. Numerical results supporting our theoretical results are provided.