Dormand-Prince and Fehlberg are two popular Runge-Kutta embedded methods for ODE integration with adaptive stepsize.
The former one estimates the error with the lower-order method and steps forward with the higher-order method; in the latter one the two roles are reversed. I do not understand the reason for this.
Why does Fehlberg method step forward with the lower-order method, instead of stepping forward with the higher-order method which is more accurate and whose result is readily available?
Whatever is the problem that I am not seeing: I expect that Dormand-Prince method is not affected by that, since it does indeed step forward with the higher-order method
The difference is not so much what step update variant is used (in RKF), but for what order the step size controller is designed.
RKF 4(5)
In RKF, you take the 4th order error $y_4-y_5=Ch^5+…$ and compare the resulting unit step error $\frac{y_5-y_4}{h}=Ch^4$ to the error tolerance $\tau$, using some scale factor involving $y_4$ and the step update $hv_4$. The resulting formula for the step size update is in its basic form $$ h_{new}=\left(\frac{\tau·h_{old}}{|y_4-y_5|}\right)^{\frac14}·h_{old}. $$
Selecting $y_5$ for the next point might get problematic if the stability region for the 4th order method is much larger than the stability region for the 5th order method. If that is the case, for a somewhat stiff ODE the step size controller will move the step size towards the boundary of the stability region of the 4th order method and oscillate about it. Then the 5th order step becomes unstable, far outside the stability region.
DoPri (4)5
In DoPri, you take the same difference $y_4-y_5=Ch^5+…$, but use the unmodified difference as proxy of the unit step error of the 5th order method. $$ h_{new}=\left(\frac{\tau}{|y_4-y_5|}\right)^{\frac15}·h_{old}. $$
This is in general wrong, the value difference is more-or-less proportional to the 5th derivative $y^{(5)}$, while the error of the 5th order step is proportional to the 6th derivative $y^{(6)}$ of the solution. So what you would actually need is more like the derivative of the leading coefficient of the 4th order error estimate, multiplied by $h^5$.
This idea of the leading error coefficient being a derivative is however only largely true for methods of order up to 4 where the order is equal to the number of stages. For higher order methods the leading coefficients of the error are polynomial expressions in the derivatives of the ODE function with different coefficients than those in the derivatives of the solution.
This intentionally wrong use of the 4th order error gets compensated in the DoPri methods by a careful construction of the method coefficients so that at least the control at the boundary of the stability region gives results suitable for the 5th order method. Another measure are the diverse filtering approaches, using the history of step sizes and error (guidance) estimates or leading coefficients of the error to get an improved, smoothed optimal step size.