I am using Boost ODEINT C++ to solve a simple differential equation with two different steppers: "runge_kutta_cash_karp54" and "runge_kutta_fehlberg78". According to the Boost website, "runge_kutta_cash_karp54" has a 5th-order error estimation and "runge_kutta_fehlberg78" has an 8th-order error estimation, which means "runge_kutta_fehlberg78" should be more accurate. However, when I solve a differential equation like "dy/dx=pxsin(x)", "runge_kutta_cash_karp54" produces more accurate results than "runge_kutta_fehlberg78" which I have compared with Mathematica result. On the other hand, for a differential equation like "dy/dx=1/(p-x)" with a singularity at "x=p", "runge_kutta_cash_karp54" throws an error in the limit of [-5,5] for p=1, which is expected, but "runge_kutta_fehlberg78" does not throw an error. I am very confused about this behaviour and would appreciate any suggestions for the detailed documentation of each individual stepper.

I am attaching the code.