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). 
Numerically solving a system of stiff nonlinear PDEs
pde;numerical analysis;nonlinear equations;nonlinear programming;stiffness