Dear (more advanced) users of Mathematica,

I'm still a beginner and currently trying to solve a system of delayed partial differential equations. Two of them describe the evolution of of two optical fields ($s_\text{P}$ and $s_\text{St}$) along a waveguide (position z):

(1) $\frac{\partial s_\text{P}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R\left(z,t\right) s_\text{St}$, and (2) $\frac{\partial s_\text{St}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R^*\left(z,t\right) s_\text{P}\left(z,t\right)$.

These equation are coupled by the amplitude $R\left(z,t\right)$ of vibrations driven by the optical fields, which itself depends on time:

(3) $\frac{\partial R\left(z,t\right)}{\partial t} + \frac{\Gamma}{2} R\left(z,t\right) = i \kappa_\text{R} s_\text{P}\left(z,t\right) s_\text{St}^* \left(z,t\right)$.

I can obtain the steady-state solution ($\partial R/\partial t =0$) very easily by putting (3) into (1),(2) and applying NDSolve and the boundary conditions ($s_0\left(0,t\right) = \sqrt{0.9}$ and $s_\text{-1}(0,t)=\sqrt{0.1}$) to these two equations. But what I'd be interested in is the build-up of R over time. If I simply integrate (3) with $R(z,0) = 0$, I end up with an "delayed PDE" error.

Does anyone have an idea how I could proceed?

Thanks for your help in advance!