eqns = {D[g[z, t], z] == 
   0.1 (\[Alpha] (1 - 1/(sqrt[1 + a^2/2]) + 
          a^2*(1 + g[z, t]/f[z, t]^2)^(-3/2)) - 
       1/(k^2*r0^2*f[z, t]^2))*(D[g[z, t], 
        t] - ((2*g[z, t]/f[z, t])*D[f[z, t], t])) - ((4*g[z, t]*
         Rd)/(t0*k^2*r0^2*f[z, t]^3))*D[f[z, t], t],
  D[f[z, t], z, z] == ( 1/(f[z, t]^3)) - 
    0.5*(\[Omega]p\[Omega])^2*(k^2*
       r0^2/\[Epsilon]0)*(g[z, 
        t]/((1 + g[z, t]/2 f[z, t]^2)^(3/2)*f[z, t]^3))}
bcs = {g[0, t] == Exp[-t^2], f[0, t] == 1, 
   Derivative[1, 0][f][0, t] == 0};
ics = {g[z, 0] == 1, f[z, 0] == 1 };

sol = NDSolveValue[{eqns, bcs, ics}, {g, f}, {z, 0, 5}, {t, -3, 3},
  Method -> {"IndexReduction" -> Automatic,
    "EquationSimplification" -> "Residual",
    "PDEDiscretization" -> {"MethodOfLines",
      "SpatialDiscretization" -> {"TensorProductGrid",
        "MinPoints" -> 200, "MaxPoints" -> 200}}}]

I have used the technique of index reduction method. I followed one example from documentation in which there is system of pde. One equation have second order spatial derivative while other have 1st order time and spatial derivation both. My problem is similar to that documentation problem. However i am not able to get any answer till now. Please find full program in the attachment. Please someone help.

Attachments: