differential equations - Numerical solution to PDE seems to contradict initial values

I am trying to solve this PDE system:

EPS = NDSolveValue[{
    D[e[z, t], z] == 0,
    D[p[z, t], t] == -p[z, t] + I s[z, t],
    D[s[z, t], t] == I p[z, t],
    e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0

    }, {e, p, s}, {z, 0, 1}, {t, 0, 1}];

Then I seek values for e and p at different z and t=0:

{EPS[[1]][0, 0], EPS[[2]][0, 0]}

{EPS[[1]][1, 0], EPS[[2]][1, 0]}

The result contradicts the condition p[z, 0] == 0:

{1, I}
{1, I}

Furthermore, if I ONLY change the order of the equations

EPS = NDSolveValue[{
    D[p[z, t], t] == -p[z, t] + I s[z, t],
    D[s[z, t], t] == I p[z, t],
    D[e[z, t], z] == 0,
    e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0
    }, {e, p, s}, {z, 0, 1}, {t, 0, 1}];

{EPS[[1]][0, 0], EPS[[2]][0, 0]}


{EPS[[1]][1, 0], EPS[[2]][1, 0]}

The result changes, and is wrong again

{1, 0.999909 I}
{1.76683*10^10, - 1.09719*10^6 I}

Does anyone have an idea what is happening?

UPDATE: The original problem arised from this fully coupled equation:

NDSolveValue[{

D[e[z, t], z] == I p[z, t],
D[p[z, t], t] == -p[z, t] + I e[z, t] + I g[t] s[z, t],
D[s[z, t], t] == I Conjugate[g[t]] p[z, t],
s[z, 0] == p[z, 0] == 0, e[0, t] == f[t]
}, {e, p, s}, {z, 0, 1}, {t, 0, 1}]

Here g[t] and f[t] are simple functions, such as const, Gaussian, etc.

Answer

The first thing I tried was to check if this is actually time integrated and not treated as a pure spatial PDE. Just calling a variable t does not tell NDSolve that something is to be considered time dependent.

EPS = NDSolve[{

    D[e[z, t], z] == 0,
    D[p[z, t], t] == -p[z, t] + I s[z, t],
    D[s[z, t], t] == I p[z, t],
    e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0, 
    1}, {t, 0, 1}, Method -> "MethodOfLines"][[1]]

NDSolve::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.

OK, so the PDE is treated as a pure spatial problem and it used the FEM for that:

EPS = NDSolve[{

     D[e[z, t], z] == 0,
     D[p[z, t], t] == -p[z, t] + I s[z, t],
     D[s[z, t], t] == I p[z, t],
     e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0, 
     1}, {t, 0, 1}][[1]];
(e /. EPS)["ElementMesh"]

NDSolve`FEM`ElementMesh[{{0., 1.}, {0., 
   1.}}, {NDSolve`FEM`QuadElement["<" 400 ">"]}]

Now, we look at the equation once more and when we eliminate the constant derivative condition we get:

EPS = NDSolve[{
    (*D[e[z,t],z]\[Equal]0,*)
    D[p[z, t], t] == -p[z, t] + I s[z, t],
    D[s[z, t], t] == I p[z, t],
    e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0, 
    1}, {t, 0, 1}][[1]]
{e -> Function[{z, t}, 1], 
 p -> InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>], 
 s -> InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>]}


{(e /. EPS)[0, 0], (p /. EPS)[0, 0]}
{(e /. EPS)[1, 0], (p /. EPS)[1, 0]}
{1, 0.}
{1, 0.}

This time a time integration did happen and for the spatial discretization you can use either FEM or TGP by givin the option (if you want): Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid"}} or "FiniteElement".

In the pure spatial problem a slightly different PDE is solved. Now, the reordering is an unfortunate issue that is documented and there is currently no way to reliably autodetect that and warn about it. For an explanation and workarounds see this tutorial in the section about Ordering of Dependent Variable Names.