High order numerical method for a subdiffusion problem

We consider a subdiffusive fractional differential problem characterized by an equation that incorporates a time Riemann-Liouville fractional derivative of order $1-\alpha$, $\alpha \in (0,1)$, on its right-hand side, while the diffusive coefficient is allowed to vary with both space and time. An high order numerical method for the subdiffusion problem is derived based on the fractional splines of degree $\beta \in (1,2]$. The main purpose of this work is to apply fractional splines for approximating the fractional integral in the definition of the Riemann-Liouville fractional derivative, and hence explain how to solve the subdiffusion problem using this approach. It is discussed how the rate of convergence of the numerical method depends on the solution, the degree of the spline and the order of the fractional derivative.