Abstract We develop two kinds of numerical schemes to efficiently solve the subdiffusion equation, which is used to describe anomalous subdiffusive transport processes. The time fractional derivative is first discretized by L1-approximation and the Grünwald–Letnikov approximation, respectively. Then we use the orthogonal spline collocation method to approximate the two semi-discretized subdiffusion equations. The stability and convergence of time semi-discretization and full discretization schemes are both established strictly for the two schemes. Both of them are unconditionally stable. Numerically the convergent orders in space (including the solution and its first derivative) are four for the Hermite cubic spline approximation, and theoretically we get that at least the solution itself has a fourth order convergent rate. Extensive numerical results are presented to show the convergent order and robustness of the numerical schemes.