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_1beta 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(icdot k_zcdotDelta 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(-icdot k_zcdotDelta 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.5Delta 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)=-hcdot f^3$$
on WolframAlpha I get
$$f(t)=pmfrac{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_zhat{E}(omega, z)=ik_zhat{E}(omega, z)-frac{omega}{2varepsilon_0 c^2k_z}hat{J}(omega, z)$$
with $J$ given by the non-linear $K$-photon-absorption ($K=2$):
$$begin{split}J&=frac{2varepsilon_0n_0csigma_2I^2hbaromega_0varrho_{nt}}{I}E\
&=varepsilon_0^2n_0^2c^2sigma_2hbaromega_0E^3end{split}$$
with $sigma_2$ the absorption cross section which can be reduced to the equation given above:
$$J=c_1beta 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 $Jgeq 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.