I am trying to solve a system of equations and have a question regarding the validity of my approach when implementing a fifth-order Cash-Karp Runge-Kutta (CKRK) embedded method with the method of lines. To give the questions some context, let me state the problem I am attempting to solve:$$ \frac {\partial E}{\partial z} = - \frac {1}{c^2}\frac {\partial E}{\partial t} - \frac{1}{k} \frac {\partial ^2 E}{\partial z^2} - \frac{1}{kc^2} \frac {\partial^2 E}{\partial t^2} + iP \tag{1}$$ $$\\\frac {\partial P}{\partial t} = iNE^* \tag{2}\\$$ $$\frac {\partial N}{\partial t} = iPE \tag{3}$$ $$ E(z=0) = \frac{\partial E}{\partial z}(z=0) = E(t=0) = \frac{\partial E}{\partial t}(t=0) = 0,\\ P(t=0) = P_0e^{z/c}, N(t=0) = N_0e^{z/c} $$where $ c = 3 \times 10^8, k = 1000, P_0 $ and $N_0$ are constants$ i=\sqrt{-1}$$;  0 \leq t \leq 1000, 0 \leq z \leq 1000 $I am implementing CKRK on the above, and even though the first spatial derivative of $E$ depends on the second spatial derivative of $E$, the numerical method appears to work when solving (1)-(3) when I use the scheme of approximating the time and spatial derivatives of $E$ on the right hand side of (1) by a backward difference approximation (I am using an accuracy of 5).To switch (1) above to a system of first order spatial derivatives in $z$, I could make the substitution:$$ U = \frac {\partial E}{\partial z} $$and solve the following equations for $E$ instead:$$ U = \frac {\partial E}{\partial z} \tag{4}\\$$ $$\frac {\partial U}{\partial z} = - \frac {k}{c^2}\frac {\partial E}{\partial t} - kU -  \frac{1}{c^2}\frac {\partial^2 E}{\partial t^2} + kP \tag{5}$$But when testing these same initial/boundary conditions using the same numerical method on the coupled equations (2) - (5), the code takes too long to finish (the step sizes required become extremely small). I believe it is due to the fact that the coefficients on the right hand side of (5) are very large and cause stability issues. I have tried to rescale the values for $z,t,P,N,$ and $E$, but doing so causes one of the other coupled equations to become unstable or has no effect (e.g. scaling $z$ does nothing to the value $ E = U\Delta z$ since both $U$ and $\Delta z$ would scale reciprocally and cancel any effect). It is due to similar reasons I am solving $E$ in the $z$-direction as opposed to doing the substitution $U = \frac {\partial E}{\partial t}$ and solving it in $t$ which is the standard method of lines approach (when I tried this method, the $\Delta t$ given by CKRK becomes very small).So ultimately, instead of using equations (2) - (5), I was wondering if applying CKRK to (1) - (3) is still a valid approach where I approximate the derivatives of $E$ on the right-side of (1) by backward finite differences? It seems very odd to apply CKRK to a first order spatial derivative that depends on an approximation of the second order spatial derivative, but is this wrong? (I would be using stored intermediate values of $E$ to ensure the backward finite difference approximations are also following the Runge-Kutta method.) 
Using backward difference approximations for higher order derivatives
pde;numerical analysis;runge kutta;method of lines