“…These models exhibit multi-scaling time behaviour, which makes them suitable for the description of different diffusive regimes and characteristic crossover dynamics in complex systems [21,25,30]. They are always formulated in the integral form, including the Riemann-Liouville fractional integral (I β t w)(t) := In capturing the multi-scale behaviors in many of integro-differential equations, such as the timefractional phase field models [8,9,14,15,[31][32][33]35] and nonlinear fractional wave models [1,2,[4][5][6]19,20,26], adaptive time-stepping strategies, namely, small time steps are utilized when the solution varies rapidly and large time steps are employed otherwise, are practically useful [3,15,18,20,22,23,[26][27][28][29]. It requires practically and theoretically reliable (stable and convergent) time-stepping methods on general setting of time step-size variations [8,9,14,15,18].…”