In this paper we aim at developing highly accurate and stable method in temporal direction for time-fractional diffusion equations with initial data $$v\in L^2(\varOmega )$$ . To this end we begin with a kind of (time-)fractional ODEs, and a hybrid multi-domain Petrov–Galerkin spectral method is proposed. To match the singularity at the beginning of time, we use geometrically graded mesh together with fractional power Jacobi functions as basis on the first interval. Jacobi polynomials are then chosen to approximate the solution in temporal direction on the intervals hereafter. The algorithm is motivated by the discovery that the solution in temporal direction possesses high regularity in proper weighted Sobolev space on a special mesh in piecewise sense. Combining standard finite element method, we extend it to solve time-fractional diffusion equations with initial data $$v\in L^2(\varOmega )$$ , in which the mesh in temporal direction is determined by spatial mesh size and finial time T. Numerical tests show that the scheme is stable and converges exponentially in time.