# RK4-method starts oscillating above certain input parameters

I am trying to solve an equation of the following type $$\partial_zE(z)=-c_0J$$ with $$J=c_1\beta E^3(z)$$ using the boost::odeint-framework and a fixed time stepper, with $$c_0$$, $$c_1$$ and $$\beta$$ constant values.
The right hand side of the function is calculated by the following procedure (depending on the step size $$\Delta z$$):

• Multiplication of $$E$$ with a vector $$\exp(i\cdot k_z\cdot\Delta z)$$ to move the field from $$z=0$$ to $$z=\Delta z$$
• Calculation of $$J$$, with the equation given above
• Multiplication of $$J$$ with a vector $$\exp(-i\cdot k_z\cdot\Delta z)$$ to move the result back from $$z=\Delta z$$ to $$z=0$$
• Return the value to the solver

This means that for a common RK4-approach I have to call that function four times, once for $$z=0$$, twice for $$z=0.5\Delta z$$ and once for $$z=\Delta z$$. This approach works fine for $$\max(E)\ll E_{lim}$$, but for $$\max(E)\approx E_{lim}$$ I have to reduce my step size $$\Delta z$$ significantly, else my resulting values for $$E$$ will explode within a few steps. $$E_{lim}$$ can be estimated by using $$J_{lim}\approx 10^{13}$$ (which might depend on the values for $$c_0$$ and $$c_1$$). If $$J$$ stays below that threshold, even large steps are not an issue, but the closer I get to that limit, the smaller I have to choose my steps (which in turn increase my calculation time).
Why does my equation behave like that? When calculating the solution for a similar type equation $$f'(t)=-h\cdot f^3$$ on WolframAlpha I get $$f(t)=\pm\frac{1}{\sqrt{c_0+2ht}}$$

To address the comment from @MaximUmansky:
I am trying to solve a reduced version of the UPPE, written as $$\partial_z\hat{E}(\omega, z)=ik_z\hat{E}(\omega, z)-\frac{\omega}{2\varepsilon_0 c^2k_z}\hat{J}(\omega, z)$$ with $$J$$ given by the non-linear $$K$$-photon-absorption ($$K=2$$): $$\begin{split}J&=\frac{2\varepsilon_0n_0c\sigma_2I^2\hbar\omega_0\varrho_{nt}}{I}E\\ &=\varepsilon_0^2n_0^2c^2\sigma_2\hbar\omega_0E^3\end{split}$$ with $$\sigma_2$$ the absorption cross section which can be reduced to the equation given above: $$J=c_1\beta E^3(z)$$ For solving the equation I can split it into the linear and the non-linear part, and start by solving the non-linear part using a RK4-solver, and applying the linear part afterwards. Implementation of the right-hand side function is similar as described above, except by the addition of an FFT for transferring $$\hat{E}\rightarrow E$$ and back, before and after the calculation of $$J$$.

Still, as soon as $$J$$ reaches a value $$J\geq J_{lim}$$ (which occurs for either large values of $$E$$ or larger values for $$\sigma_2$$, i.e. a large absorption), the problems described above occur, and the closer $$J$$ gets to $$J_{lim}$$, the smaller I have to choose my step size, to prevent $$J$$ becoming larger than $$J_{lim}$$, which is a behavior I still do not understand.

• Is this from electrodynamics in continuous media? It would help to know the physical meaning of those equations, that may help understand the problem you see in the calculations. – Maxim Umansky Jul 16 at 15:42
• @MaximUmansky: I extended my question, adding additional information. If those are not sufficient, please let me know. – arc_lupus Jul 17 at 9:05
• Have you tried to employ a more stable method than RK4? Most of the stable method are implicit, but there exist explicit alternative. You can have a look at Runge-Kutta-Chebyshev (RKC) methods, or ROCK. – G. Gare Jul 17 at 9:37
• @G.Gare: I'd like to avoid implicit methods, usually my matrix for $E$ has $>1e6$ entries, but will take a look at RKC and ROCK methods. – arc_lupus Jul 17 at 10:05