Over the past twenty years, the mathematical modeling of Fractional Differential Equations (FDEs) has attracted growing attention and expanded significantly across various scientific fields, including bioengineering, mathematical biology, physics, chemistry and computational medicine. In this context, finding analytical solutions of FDEs is often more challenging than for classical ordinary differential equations, while the derivation of accurate and reliable numerical methods suffers from the possible non-smoothness of the solution and/or the vector field at the starting time, not to mention that the efficient treatment of the persistent memory term can make long-time simulations computationally demanding, due to the non-locality of the operator. To mitigate the aforementioned issues, the class of Runge–Kutta-type methods, named Fractional HBVMs (FHBVMs) is presented, covering its design, development and analysis. In particular, a novel extension is proposed, allowing for a mixed graded/uniform stepsize strategy that gives rise to an updated version of the pre-existing fhbvm Matlab code. As confirmed by the numerical tests, the novel approach is mainly tailored to problems with a non-smooth vector field at the starting time, whose solution is both non-smooth at the origin and oscillatory, to concurrently gain accuracy in reproducing the initial non-smooth behaviour of the solution, while maintaining efficiency over extended time periods.