Using NDSolve and PieceWise for boundary conditions for coupled PDEs

Question

I've decided to try and use PieceWise in my boundary conditions for my 4 coupled PDEs, but it doesn't seem to work out nicely. I'm receiving the error:

NDSolve::overdet: There are fewer dependent variables, {IPB[z,t],IPF[z,t],IS1B[z,t],IS1F[z,t]}, than equations, so the system is overdetermined.

Paramters:

ROCP = 1;
RICP = 0;
ROCS1 = 0.5;
RICS1 = 1;
L = 5 10^-2;
n = 1.556;
c = 3 10^10;
lR = 3;
R1 = 0.99;

gP = 20;
gS1 = 20;
sigmaSP = 10^-13;

toffset = 5 10^-9;
twidth = 10 10^-9;
Energy = 1.3 10^-13;
rad = 0.015;
PumpInt = ((Energy)/(twidth ))/(Pi rad^2) ;
PumpPeak = Exp[0.5 ((-toffset)/(twidth))^2] PumpInt;

lc = n lR;
alphaP = L/(2 lc);
alphaS1 = (L - Log[R1])/(2 lc);

tmax = 20 10^-9;

The 4 PDEs I'm solving are:

PDE = {D[IPF[z, t], z] == -(n/c) D[IPF[z, t], t] - 
     gP IPF[z, t]*(IS1F[z, t] + IS1B[z, t]) - alphaP IPF[z, t],
   -D[IPB[z, t], z] == -(n/c) D[IPB[z, t], t] - 
     gP IPB[z, t] (IS1F[z, t] + IS1B[z, t]) - alphaP IPB[z, t],
   D[IS1F[z, t], z] == -(n/c) D[IS1F[z, t], t] + 
     gS1 IS1F[z, t] (IPF[z, t] + IPB[z, t]) - alphaS1 IS1F[z, t] + 
     sigmaSP (IPF[z, t] + IPB[z, t]),
   -D[IS1B[z, t], z] == -(n/c) D[IS1B[z, t], t] + 
     gS1 IS1B[z, t] (IPF[z, t] + IPB[z, t]) - alphaS1 IS1B[z, t] + 
     sigmaSP (IPF[z, t] + IPB[z, t])};

And the Boundary conditions are:

BC = {IPF[z, t] == 
    PieceWise[{{PumpPeak*Exp[-0.5 ((t - toffset)/twidth)^2], 
       z == 0 && t > 0}, {0, 0 <= z <= lR && t == 0}}],
   IPB[z, 0] == 0, IPB[lR, t] == IPF[lR, t] ROCP,
   IS1F[0, t] == IS1B[0, t] RICS1, IS1B[lR, t] == IS1F[lR, t] ROCS1, 
   IS1F[z, 0] ==  0, IS1B[z, 0] ==  0};

And the NDSolve:

solInt = NDSolve[{PDE, BC}, {IPF, IPB, IS1F, IS1B}, {z, 0, lR}, {t, 0,
     tmax}, MaxStepFraction -> {1/750, 1/10}, 
   Method -> {"BDF", "MaxDifferenceOrder" -> 5}, 
   InterpolationOrder -> All];

The main idea behind the boundary conditions is that IPF[z,t] = 0 at time t = 0 and for any z in the range 0 - lR, but for t > 0 and z = 0, it will have the form of the Gaussian equation.

Is it possible to impose these initial conditions?

Side note: Sorry if my post are similar, I didn't know if I should piggyback on my last question or if this was something I should ask from scratch.

Cheers!

Follow up: I would just like to try and extend this question since I feel this may be related. After modifying the parameters and boundary and initial conditions I now have:

ROCP = 1;
RICP = 0;
ROCS1 = 0.5;
RICS1 = 1;
L = 5 10^-2;
n = 1.556;
c = 30;
lR = 3;
R1 = 0.99;

gP = 20 ;
gS1 = 20 ;
sigmaSP = 10^-13;

toffset = 10;
twidth = 2.5;
Energy = 1.3 10^-13;
rad = 0.015;
PumpInt = ((Energy)/(twidth 10^-9))/(Pi rad^2) 
PumpPeak = Exp[0.5 ((-toffset)/(twidth))^2] PumpInt

lc = n lR;
alphaP = L/(2 lc);
alphaS1 = (L - Log[R1])/(2 lc);

tmax = 20;
ppR = 100;

And conditions:

    BC = {IPF[0, t] == PumpPeak*Exp[-0.5 ((t - toffset)/twidth)^2], 
    IPB[lR, t] == IPF[lR, t] ROCP, IS1F[0, t] == IS1B[0, t] RICS1, 
    IS1B[lR, t] == IS1F[lR, t] ROCS1,
    IPB[z, 0] == 0, IPF[z, 0] == 0, IS1F[z, 0] ==  0, 
    IS1B[z, 0] ==  0}

But when I go to solve:

solInt = 
  NDSolve[{PDE,BC},
   {IPF, IPB, IS1F, IS1B}, {z, 0, lR}, {t, 0, tmax}, 
   MaxSteps -> Infinity, PrecisionGoal -> 2, AccuracyGoal -> 50, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> ppR, "MaxPoints" -> ppR, 
       "DifferenceOrder" -> 2}, 
     Method -> {"Adams", "MaxDifferenceOrder" -> 1}}];

I get the errors:

NDSolve::ndsz: At t == 0.6079745549156106`, step size is effectively zero; singularity or stiff system suspected.

NDSolve::eerr: Warning: scaled local spatial error estimate of 30.865283175219346at t = 0.6079745549156106 in the direction of independent variable z is much greater than the prescribed error tolerance. Grid spacing with 100 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

I've narrowed it down to the term twidth*10^-9 in the PumpInt paramter as when I change it to twidth*10^9 the code runs smoothly. Unfortunately for me however, I need the factor of 10^-9. My question is how can I make the code run with this factor. Thank you for all the help so far.

Update: It could potentially be sigmaSP, since setting that term to 0 resolves the issue. Is there anyway I can include this term?