I am implementing an ODE solver using the Cash-Karp method on equations with the following form:$$ \frac {d E}{d z} = - \frac {1}{\mu_0c} \frac {d ^2 E}{d z^2} + \frac {i}{\mu_0c}E \tag{1} $$And the output appears to be working fine (the second spatial derivative being approximated by backward finite differences). Although when I convert the system to a system of first order ODEs:$$ U = \frac {\partial E}{\partial z} \tag{2}$$$$ \frac {d U}{d z} = - \mu_0c U + iE \tag{3} $$When using Cash-Karp on this new system, here is an image of the amplitude vs z output.. The amplitude displays unwanted high frequency oscillations. These oscillations should not be present and are resulting in the adaptive step becoming very slow since I am tackling larger domains and varying parameters/boundary conditions yet these oscillations remain constant at a spatial frequency of $50$ $\text{m}^{-1}$. I suspect the oscillations are due in part because the constant $\mu_0c \approx 100$ is large, yet a constraint is that I cannot rescale the variables/domain (see: full explanation of problem).I have tried adjusting the adaptive step method and varying the tolerances, but these fast oscillations appear to remain. Any suggestions on how to remove these numerical oscillations while using an explicit scheme would be of interest.
Preventing numerical oscillations with Cash-Karp method
runge kutta;numerical modelling;adaptive mesh refinement;oscillations;numerical limitations