The constants that figure in the step selection heuristics are usually obtained from numerical experiments and are deemed to be "good practice". For the method to be robust, beside this, you need to have a step rejection criterion (for example if a current LTE is larger from the previous one by a certain factor).

Thank you! Yes, there is a step rejection criterion in ode45 (if err > rtol, then the step is failed). In this case the step is changing according to the formula

absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));

As I understand, here constants 0.1 and 0.8 were obtained from experiments too, as you said.

As for the pow degree, the ratio err/tol is proportional to h^5, so pow equals 1/5, this is clear to me.

the idea is to mantain err < rtol, but we want a fast integration so we need a stepsize, h, that makes err = rtol. Or err/rtol = 1.

In the other hand, we know the error estimation is something like C*h^5 + O(h^6) where C is the main coefficient. If we reject the rest of the expresion we can say that with h, the actual stepsize, err/tol \approx C*h^5. If the stepsize were optimun, we call h_o to optimun stepsize, err/tol would be 1, so C \approx 1/h_o^5. We put this C in the previous and we obtain:

err/tol \apporx (1/h_o^5) * h^5

solving for h_o ( I will use = and not \approx):

h_o^5 = h^5 /(err/tol)

and finally,

h_0 = h / ((err/tol)^(1/5))

this is your absh = absh/temp.

The problem is that we are playing with fire, because we have obtained the optimun stepsize that would get err/tol=1. This is dangerous, because if err/tol > 1 the step will be rejected (lost CPU time!), so we have to be carefull. A way to avoid rejections is to be conservative and avoid big stepsize changes, in your case there are two controls:

1.- you do not use (err/rtol)^pow, you are using 1.25*(err/rtol)^pow

2.- the new stepsize will never be greater than 5* old stepsize ( if (temp < 0.2) absh = 5.0*absh).

I hope this will help you to understand stepsize selection.

2 Joseba Makazaga

Ok, thank you! I cleary understand how stepsize selection works and what does powers 1/5 and 5 mean, but my question was about numbers like 1.25, 0.2 and 5.0. For example, why 1.25, but not 1.35 or 2. Apparantly, they are obtained experimentally.

2 Sebastian Schöps

Tank you! There are some ideas on how this safety factor (1.25) can be calculated (for implicit Runge-Kutta, though).