I discovered a bug; your error estimation scheme (the two deltas) is incorrect! Take a look at my RKN1210 integrator how to do it; it will significantly improve your integrator's performance.
Theoretically, if you copy-paste your coefficients into my integrator and adjust the for loop limits from 1:17 to 1:8 and pow=1/12 into pow=1/8, your integrator should inherit all the additional features that I implemented :) I haven't tried this, but I think it's worth a shot.
22 Sep 2004
ode86
ode86 integrates a system of ordinary differential equations at stringent tolerances
Author: Charalampos Tsitouras
., drusi
Doesn't accept params to "evalfunc"
10 May 2004
ode86
ode86 integrates a system of ordinary differential equations at stringent tolerances
Author: Charalampos Tsitouras