In the present study, we introduce a high-order non-polynomial spline method designed for non-linear time-fractional reaction–diffusion equations with an initial singularity. The method utilizes the L2-1σ scheme on a graded mesh to approximate the Caputo fractional derivative and employs a parametric quintic spline for discretizing the spatial variable. Our approach successfully tackles the impact of the singularity. The obtained non-linear system of equations is solved using an iterative algorithm. We provide the solvability of the novel non-polynomial scheme and prove its stability utilizing the discrete energy method. Moreover, the convergence of the proposed scheme has been established using the discrete energy method in the L2 norm. It is proven that the method is convergent of order min{θα,2} in the temporal direction and 4.5 in the spatial direction, where α∈(0,1) denotes the order of the fractional derivative and the parameter θ is utilized in the construction of the graded mesh. Finally, we conduct numerical experiments to validate our theoretical findings and to illustrate how the mesh grading influences the convergence order when dealing with a non-smooth solution to the problem.