I am attempting to numerically solve the following:
\begin{align} \frac {\partial y_1}{\partial t} &= i(y_2y_3 - y_2^*y_3^*) - y_1 \tag{1}\\ \frac {\partial y_2}{\partial t} &= y_1^*y_3 - y_2 \tag{2}\\ \frac {\partial y_3}{\partial t} &= y_4 \tag{3} \\ \frac {\partial y_4}{\partial t} &= y_4+ \frac {\partial^2 y_3}{\partial z^2} + \frac {\partial y_3}{\partial z} + \frac {\partial^2 y_2}{\partial z^2} + \frac {\partial y^*_2}{\partial z} + y_1^*y_4 + y_3[i(y_2y_3 - y_2^*y_3^*) - y_1] - iy_1y_3 -y_2 \tag{4} \end{align}
The initial/boundary conditions are:
\begin{align} y_1(z,0) = 10^{-5}; y_2(z,0) = 10^4; y_3(z,0) = 0; y_4 (z,0) = 0; \end{align} \begin{align} y_1(0,t) = 10^{-5}e^{-t}; y_2(0,t) = 10^{4}e^{-t}; y_3(0,t) = 0; \frac {\partial y_3}{\partial z} (0,t) &= 0 \end{align}
With a current approach using method of lines, I implemented a central difference scheme for the spatial derivatives and Runge-Kutta (4th order) for the time derivatives. This worked but only for particular conditions within a narrow time range (after which the output blew up).
I have attempted using backward Euler instead of RK4 for each of the equations, although since it involves complex conjugates, I end up needing to split the real and imaginary parts which then yields a system of 8 equations. Splitting the real and imaginary parts is no issue itself, however, since these equations are nonlinear, I cannot isolate all the desired variables and insert them sequentially into the last equations.
Do you have any suggestions on how to tackle this problem? Any better-suited explicit methods I should try for this problem? Advice on formulating/implementing an implicit method to solve the above equations? I am using Python and have looked into existing packages, although have had no luck with a couple (e.g. FiPy).