System of nonlinearly coupled PDEs

Similar to the question I asked on the Math Stackexchange (in this question), I am interested in solving a nonlinearly coupled system of PDEs.

The system of PDEs look like this:
$$begin{align*}partial_t p_{11} &= E_p p_{13} – E_p p_{31}
partial_t p_{12} &= E_p p_{13} – E_p p_{32} – p_{12} (Delta_{31}-Delta_{32})
partial_t p_{13} &= E_p p_{11} + E_c p_{12} – E_p p_{33} + p_{13} Delta_{31}
partial_t p_{21} &= E_p p_{23} – E_c p_{31} + p_{21}(-Delta_{31} + Delta_{32})
partial_t p_{22} &= E_c p_{23} – E_c p_{32}
partial_t p_{23} &= E_p p_{21} + E_c p_{22} – E_c p_{33} + p_{23} Delta_{31} + p_{23} (-Delta_{31} + Delta_{32})
partial_t p_{31} &= – E_p p_{11} – E_c p_{21} + E_p p_{33} – p_{31} Delta_{31}
partial_t p_{32} &= – E_p p_{21} – E_c p_{22} – E_c p_{33} + p_{23} Delta_{31} + p_{23} (-Delta_{31} + Delta_{32})
partial_t p_{33} &= E_p p_{13} – E_c p_{23} + E_p p_{31} + E_c p_{32}
end{align*}$$

$$
frac{partial E_p}{partial z } + frac{1}{c} frac{partial E_p}{partial t } = i k p_{13}(t, z)
$$

Each of these p’s ($p_{11}, p_{12}, …, p_{33}$) are a function of t and z. Additionally $E_p$ is a function of t and z. And $Delta_{31, 32, 33}$ are just numerical constants.

I’m trying with the boundary conditions:
$$text{p11}(0,z)=1,
text{p12}(0,z)=text{p13}(0,z)=text{p21}(0,z)= text{p22}(0,z)=text{p23}(0,z)=text{p31}(0,z)= text{p32}(0,z)= text{p33}(0,z)=0,
E_p(t,0)=e^{-(t-t_0)^2}- e^{-(t_0)^2}
E_p(0,z)=0
$$

Is there a way that I can solve this numerically in Mathematica?

When I try the following code:

NDSolve[{
  Derivative[1, 0][p11][t, z] == -(1/2) I Ep[t, z] (p13[t, z] - p31[t, z]), 
  Derivative[1, 0][p12][t, z] == 
      1/2 I (-10 p13[t, z] + Ep[t, z] p32[t, z]), 
  Derivative[1, 0][p13][t, z] == -(1/2) p13[t, z] - 
        1/2 I (10 p12[t, z] + Ep[t, z] (p11[t, z] - p33[t, z])), 
  Derivative[1, 0][p21][t, z] == -(1/2) p21[t, z] - 
        1/2 I (Ep[t, z] p23[t, z] - 10 p31[t, z]), 
  Derivative[1, 0][p22][t, z] == -5 I (p23[t, z] - p32[t, z]), 
  Derivative[1, 0][p23][t, z] == -(1/2)
          I (Ep[t, z] p21[t, z] + 10 p22[t, z] - I p23[t, z] - 
            10 p33[t, z]), 
  Derivative[1, 0][p31][t, z] == 
      1/2 I (10 p21[t, z] + I p31[t, z] + 
            Ep[t, z] (p11[t, z] - p33[t, z])), 
  Derivative[1, 0][p32][t, z] == 
      1/2 I (Ep[t, z] p12[t, z] + 10 p22[t, z] + I p32[t, z] - 
            10 p33[t, z]), 
  Derivative[1, 0][p33][t, z] == 
      1/2 I (10 p23[t, z] + Ep[t, z] (p13[t, z] - p31[t, z]) - 
            10 p32[t, z] + 2 I p33[t, z]), 
  Derivative[0, 1][Ep][t, z] + 
    Derivative[1, 0][Ep][t, z] == 1/2 I p12[t, z], p11[0, z] == 1, 
    p12[0, z] == 0, p13[0, z] == 0, p21[0, z] == 0, p22[0, z] == 0, 
    p23[0, z] == 0, p31[0, z] == 0, p32[0, z] == 0, p33[0, z] == 0, 
    Ep[t, 0] == E^-t^2}, {ρ11[t, z], ρ12[t, z], ρ13[t, 
      z], ρ21[t, z], ρ22[t, z], ρ23[t, z], ρ31[t, 
      z], ρ32[t, z], ρ33[t, z], Ep[t, z]}, {z, 0, 1}, {t, 0, 
    1}]

Which looks like this:

Mathematica returns the error message: NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

Any ideas?

EDIT: I modified the boundary conditions for $E_p$ as suggested in the comments. (Previously I only had E_p(t,0)=e^{-(t-t_0)^2}, and now I have added an initial condition and modified the previous B.c. to make it consistent.).