I'm developing my own generic Runge-Kutta solver, and I'm currently implementing the adaptive step-size routine. I say generic because I want to be able to test different RK implementations by only passing the solver a Butcher tableau, such as the following.\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots\\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\&b_{1}^*&b_{2}^*&\dots &b_{s}^*\\\end{array}Currently, I'm only focused on adaptive-explicit methods. As far as I understand it, a non-adaptive method would only have only the first row of $b$, and an explicit method would be zero on the diagonal and upper triangle.Essentially what I'm doing within a single iteration of my solver is passing a timestep span [$t_0, t_1$], and initial guess of the step ($h$), and an initial value of $y_0$, which could be a vector, and the function $\frac{dy}{dt} = f(t, y)$.I calculate $y(t_0 + h)$ and $y^*(t_0 + h)$ using coefficients of the Butcher tableau, calculate the element-wise error between them, and then calculate a scaling factor of the step using the maximum absolute error. The scaling factor is used to determine if $h$ is too small, a good size, or too big, and is adjusted accordingly. Currently, I'm calculating the scaling factor like this (described in this youtube video):$$ s = \frac{\epsilon \cdot h}{2 \cdot (t_1 - t_0) \cdot max(|y-y^*|)}$$where $t_0$ and $t_1$ are the beginning and end of the time step of interest, and epsilon is some error parameter set previously (e.g. 1e-8).The resulting value of $h$ would be scaled using $s$:$$h = sh$$If s is larger than 2, then I assume the error between $y(t_0 + h)$ and $y^*(t_0 + h)$ is small and I can continue with the next step. If s is less than 1, then I adjust h without moving forward and try the same timestep again. If $s$ is between 1 and 2, then I move on as well, but I don't adjust h because it's good enough.What as not readily obvious to me at the start, but now becomes so as the $tspan$ increases ($t_1-t_0$ goes up), is that the acceptable h is kept really far down. This also seems to be greatly influenced by the RK method used.For example, the following system takes about 2 seconds on my laptop with 100 pts between [0, $4\pi$]:\begin{array} {rl}y_1'(t) &= y_2 \\y_2'(t) &= -y_1 \\y_1(0) &= 0 \\y_2(0) &= 0.1 \\\end{array}Decreasing the number of sampling points to 10 within the same range ([0, $4\pi$]) now takes 16 seconds because I think the value of s is too small, due to the $t_1-t_0$ term.It may be that I'm only interested in a very few number of points (or even the final value), but this is causing the solution time to increase, even for this relatively small problem.Can someone point me to a reference that goes into more detail about this scaling ($s$) parameter? Or is there another way to solve this approach this problem?
Scaling step size in adaptive runge-kutta method
ode;runge kutta