# 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. Jul 16 '20 at 15:42
• @MaximUmansky: I extended my question, adding additional information. If those are not sufficient, please let me know. Jul 17 '20 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. Jul 17 '20 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. Jul 17 '20 at 10:05