# How to make gaussian pulse move using Boundary Conditions (one spatial dimension)?

I'm trying to solve a system of normalized coupled PDE's to model Raman scattering. Essentially, a gaussian light pulse of 300ps duration enters the raman active medium at ζ = 0 which then drives the molecules and, as a result, gets scattered to a lower frequency (1st Stokes). I'm interested in how the spatial field profile for both the main pulse and stokes pulse evolves with time. I would like to reflect the fact that the pulse is longer than the 8cm crystal and only solve the coupled PDE's within a smaller subdomain. Also, I made some attempt to get the pulse to move through the medium using Dirichlet BC's at ζ = 0 but the frames that are generated show that only the peak is affected and the pulse still remains stationary. I appreciate any help on this as I've been stuck on these issues for awhile now and have looked over Mathematica documentation many times now. Here is the code:

ClearAll
(*Physical quantities and constants*)
eps0 = 8.85*10^-12;
T2 = 25*10^-12;
c = 3*10^8;
λf = 800*10^-9;
λs = 873*10^-9;
gs = 16*10^-11;
ωl = 2*π*c/λf;
ωs = 2*π*c/λs;
tp = 300*10^-12;
pdiam = .003; (*Diameter of pump in [m]]*)
w0 = pdiam/2 ;
n[λ_] := Sqrt[
2.4069 + .01992/(λ^2 - .03773) - 0.006166*λ^2]
μl = n[.8];
μs = n[.873];
κQ = I*c*μl*μs/(16*π*ωs);
κS = I*ωs/μs;
κL = I*ωl/μl;
Γl = 3.8*10^-4;
Γs = 3.8*10^-4;
vp = 3*10^8/n[.8];
vs = 3*10^8/n[.873];
ZRL = c*tp;
vpnorm = vp*tp/ZRL;
vsnorm = vs*tp/ZRL;
T = T2/tp;
(*End of Physical Quantities and Constants*)

f[ζ_] := E^(-2 Log*ζ^2)
pdes = {D[M[τ, ζ], τ] -
vsnorm*D[M[τ, ζ], ζ] -
vsnorm*(2*ζ/(1 + ζ^2))*M[τ, ζ] +
M[τ, ζ]/T == (κQ/T)*1/(1 + ζ^2)*
EL[τ, ζ]*ES[τ, ζ]\[Conjugate],
D[ES[τ, ζ], τ] + Γs*
ES[τ, ζ] == (1/3)*κS*gs*EL[τ, ζ]^2*
ZRL*vsnorm*EL[τ, ζ]*
M[τ, ζ]\[Conjugate], (vpnorm - vsnorm)*
D[EL[τ, ζ], ζ] + (vpnorm -
vsnorm)*(ζ/(1 + ζ^2))*EL[τ, ζ] +
D[EL[τ, ζ], τ] + Γl*
EL[τ, ζ] == (1/3)*κL*gs*EL[τ, ζ]^2*
ZRL*vpnorm*M[τ, ζ]*ES[τ, ζ]};
initial = {M[0, ζ] == 0,
ES[0, ζ] == 0*E^(-2 Log*(ζ/2)^2),
EL[0, ζ] == f[ζ] };
bc = DirichletCondition[{M[τ, ζ] == 0,
ES[τ, ζ] == 0,
EL[τ, ζ] ==
E^(-2 Log*(ζ - vpnorm*τ)^2)}, ζ == 0 ];
{M, ES, EL} =
NDSolveValue[{pdes, initial, bc}, {M, ES, EL}, {τ, 0,
2.6*10^-10/tp}, {ζ, 0, .08/ZRL}];
framesTEQ =
Table[Plot[{EL[τ, ζ]*EL[τ, ζ]\[Conjugate],
ES[τ, ζ]*ES[τ, ζ]\[Conjugate]}, {ζ,
0, .08/ZRL}, PlotRange -> {0, 1}], {τ,
0, (2.6*10^-10)/tp, (2.6*10^-10)/(20*tp)}];
Manipulate[framesTEQ[[i]], {{i, 3, "time"}, 1, Length[framesTEQ], 2},
SaveDefinitions -> True]