A check that the problem is solvable in mathematica or not; But is numerically solvable;

ClearAll["Global`*"];eq=DSolve[{T*d[B[z,n],{z,2}]-k*DiracDelta[z]*B[z,n]==d[B[z,n],n],B[z,0]==Exp[-(x+a)^2/(4*s^2)]},B[z,n],z,n]
e=0.0001;
T=0.5;
k=1;
a=1;
s=10;
sol=Block[{ϵ=$MachineEpsilon},NdSolve[{T*d[B[z,n],z,z]==d[B[z,n],n],B[z,0]==Exp[-(z+a)^2/(4*s^2)],B[0,n]==0,Derivative[1,0][B][0,n]==1},B[z,n],{z,-1000,1000},{n,0,100}]]

1-Erf[x]===Erfc[x]

Plot[{1-Erf[3x/(1+4)],Erfc[3x/(1+4)]},{x,-10,10}]

Doing an errorless numerical calculation

ClearAll["Global`*"];
ϵ0=k0 /(2d);
Clear[t];
ϵ2=k0/2/d;
Q[t_]=(1+2UnitStep[-xc-x0-2ϵ2(σ^2+d t)]Exp[-ϵ2^2(σ^2+d t)]-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))]);
Qnc[t_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)](Erfc[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))])
x0=2;
k=0.1;
β=10;
σ=0.5/Sqrt[β];
h=-.4;
y=-1.5;
ys=0;
k0=3;
d=.4;
xc=1;
t1=Table[(Erfc[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))]),{t,0,50}];
t2=Table[Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)],{t,0,50}];
t3=Table[Q[t],{t,0,50}];
t1 t2+t3
p1=ListLinePlot[%,PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
p2=Plot[Qnc[t],{t,0,50},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
Show[p1,p2]

{1,2}.{2,3}

Q[t]+Q2[t]

p2=Plot[Qnc[t],{t,0,50},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}];
Table[{t,SetPrecision[Qnc[t],32]},{t,0,10,0.1}];
p3=ListLinePlot[%]
Show[p1,p2,p3]

Quantifying the error:

ClearAll["Global`*"];
ϵ0=k0 /(2d);
Clear[t];
ϵ2=k0/2/d;
Q[t_]=(1+2UnitStep[-xc-x0-2ϵ2(σ^2+d t)]Exp[-ϵ2^2(σ^2+d t)]-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))]+(1(*Sign[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))]*)-Erf[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)]);
Qnc[t_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)](Erfc[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))])
x0=2;
k=0.1;
β=10;
σ=0.5/Sqrt[β];
h=-.4;
y=-1.5;
ys=0;
k0=3;
d=.4;
xc=1;
p1=Plot[Qnc[t],{t,0,50},PlotLegends->Placed[{Style["b",13]},{Scaled[{0,0.5}],{y,h}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->Gray]
p2=Plot[Q[t],{t,0,50},PlotLegends->Placed[{Style["a",13]},{Scaled[{0,0.5}],{y,h}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Black,Dashed}]
Show[p1,p2]

ClearAll["Global`*"];
P[x_,t_]=(((Exp[((-(x+x0)^2)/(4(σ^2+d t)))])/(Sqrt[(σ^2+d t)4π]))+((ϵ0)/2)*Exp[ϵ0(x*Sign[x-xc]+2xc*UnitStep[xc-x]+x0)+ϵ0^2(σ^2+d t)]*(-Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]+Erf[(x*Sign[x-xc]+2xc*UnitStep[xc-x]+x0+2ϵ0(σ^2+d t))/(2Sqrt[(σ^2+d t)])]));
ϵ0=k0 /(2d);
Clear[t];
ϵ2=k0/2/d;
Q[t_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)](Sign[((xc+x0)/(2(σ^2+d t)))+ϵ2]-Erf[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])+2UnitStep[-xc-x0-2ϵ2(σ^2+d t)]Exp[-ϵ2^2(σ^2+d t)]-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))])
Qnc[t_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2(xc+x0)](Erfc[((xc+x0+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])-Erfc[((xc+x0)/(2Sqrt[(σ^2+d t)]))])
Pexp[x_,t_]=P[x,t];
Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
x0=2;
k=0.1;
β=10;
σ=0.5/Sqrt[β];
h=-.4;
y=-1.5;
ys=0;
k0=3;
d=.4;
xc=1;
ab1=Plot[{Pexp[x,0],Pexp[x,1],Pexp[x,2],Pexp[x,5],Pexp[x,10],Pexp[x,20],Pexp[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotTheme->"Scientific",PlotLabel->{Style["Probability distribution",14]},FrameLabel->{Style["x",14],Style["P(x,t)",14]}]
Clear[t]
p1=Plot[Q[t],{t,0,50},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}];
p2=Plot[Qnc[t],{t,0,50},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}];
Table[{t,SetPrecision[Qnc[t],32]},{t,0,10,0.1}];
p3=ListLinePlot[%]
Show[p1,p2,p3]
S[z_]=Sqrt[1000/π]Exp[-1000* (z)^2];
dThick=20;
dTime=50;
eqSimplee=D[Pe[x,t],t]==d*D[Pe[x,t],{x,2}]-k0*S[x-xc]*Pe[x,t];
eqInie=Pe[x,0]==√(β/π)*Exp[-β*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqInie,eqBlefte,eqBleftg},{Pe},{x,-dThick,dThick},{t,0,dTime}];
ab2=Plot[{Evaluate[Pe[x,0]/.res],Evaluate[Pe[x,1]/.res],Evaluate[Pe[x,2]/.res],Evaluate[Pe[x,5]/.res],Evaluate[Pe[x,10]/.res],Evaluate[Pe[x,20]/.res],Evaluate[Pe[x,50]/.res]},{x,-20,20},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->Full,PlotTheme->{"Scientific"},PlotLabel->{Style["Probability distribution",14]},FrameLabel->{Style["x",14],Style["P(x,t)",14]}]
Show[ab1,ab2]
t1=Flatten[Table[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick},MaxRecursion->20],{t,0,50,5}]]
t2=Table[{t,t1[[t/5+1]]},{t,0,50,5}]
p4=ListPlot[%,PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Dashed}];
Show[p2,p3]

Show[p2,p4]

Clear[xc];
xc=-2;
Plot[Sign[((xc+x0)/(2(σ^2+d t)))+ϵ2],{t,0,50}]

x0 = 10;β = 10;σ = 0.5/Sqrt[β];k0 = 2;d = .1;xc = 1;

Some important integral involved;

ClearAll["Global`*"];(-I*a)/(2*Pi)*Integrate[Exp[-(σ^2+d*t)u^2]/(u+I*(a+b)),{u,-Infinity,Infinity},Assumptions->Element[σ,Reals]&&σ^2+d*t>0&&Element[a,Reals]&&Element[b,Reals]&&(a+b)>0]
(-I*a)/(2*Pi)*Integrate[Exp[-(σ^2+d*t)u^2]/(u+I*(a+b)),{u,-Infinity,Infinity},Assumptions->Element[σ,Reals]&&σ^2+d*t>0&&Element[a,Reals]&&Element[b,Reals]&&(a+b)<0]

and transform calculations to compare with arxiv: 1408.3288
https://arxiv.org/pdf/1408.3288.pdf;

ClearAll["Global`*"];
P[x_,t_]=(E^(-(1/4) (((x+x0))^2) (d t+(σ^2))))-(1/4) (E^(((ϵ0)^2) (d t+(σ^2))-ϵ0 (x0+Abs[x]))) ϵ0 (1-Erf[((x0+2 ϵ0 (d t+(σ^2))+Abs[x])/(2 (Sqrt[d t+(σ^2)])))] Sign[x0-2 ϵ0 (d t+(σ^2))+Abs[x]])
Integrate[P[x,t]Exp[-s*t]/(2*Sqrt[Pi*((p^2)+d*t)]),{t,0,Infinity},Assumptions->Element[s,Reals]&&d>0&&a<0&&Element[p,Reals]]
LaplaceTransform[P[x,t],t,s]
Pe[x_,s_]=Exp[Sqrt[s/A]*x]+Erf[x0+2*p^2Sqrt[s/A]/p]
InverseLaplaceTransform[Pe[x,s],s,t]

ClearAll["Global`*"];
P[x_,t_]=(((Exp[((-(x+x0)^2)/(4(σ^2+d t)))])/(Sqrt[(σ^2+d t)4π]))+((ϵ0)/2)*Exp[ϵ0(x*Sign[x-xc]+2xc*UnitStep[xc-x]+x0)+ϵ0^2(σ^2+d t)]*(-Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]+Erf[(x*Sign[x-xc]+2xc*UnitStep[xc-x]+x0+2ϵ0(σ^2+d t))/(2Sqrt[(σ^2+d t)])]));
ϵ0=k0 /(2d);
Clear[t];
ϵ2=k0/2/d;
Q[t_,x0_,xc_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2*RealAbs[xc+x0]](Erfc[((RealAbs[xc+x0]+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])+2UnitStep[-RealAbs[xc+x0]-2ϵ2(σ^2+d t)]Exp[-ϵ2^2(σ^2+d t)]-Erfc[((RealAbs[xc+x0])/(2Sqrt[(σ^2+d t)]))]);
Q2[t_,x0_,xc_]=(1+Exp[ϵ2^2(σ^2+d t)+ϵ2*RealAbs[xc+x0]](Sign[((RealAbs[xc+x0])/(2(σ^2+d t)))+ϵ2]-Erf[((RealAbs[xc+x0]+2ϵ2(σ^2+d t))/(2Sqrt[(σ^2+d t)]))])+2UnitStep[-RealAbs[xc+x0]-2ϵ2(σ^2+d t)]Exp[-ϵ2^2(σ^2+d t)]-Erfc[((RealAbs[xc+x0])/(2Sqrt[(σ^2+d t)]))]);

Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
k=0.1;
β=10;
σ=0/Sqrt[β];
h=-.4;
y=-1.5;
ys=0;
k0=1;
d=1;
p1=Plot[{Q[t,10,0]},{t,0,50},PlotLegends->Placed[{Style["((Subscript[x, 0]),(Subscript[x, c]))=(0,1)",13]},{Scaled[{0,0.5}],{y,h}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},PlotStyle->Gray,FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
p2=Plot[{Q[t,-10,0]},{t,0,50},PlotLegends->Placed[{Style["((Subscript[x, 0]),(Subscript[x, c]))=(0,-1)",13]},{Scaled[{0,0.5}],{y,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},PlotStyle->Dashed,FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
p3=Plot[{Q[t,3,-13]},{t,0,50},PlotLegends->Placed[{Style["((Subscript[x, 0]),(Subscript[x, c]))=(2,-3)",13]},{Scaled[{0,0.5}],{y,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},PlotStyle->{Black,Dotted},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
Show[p1,p2,p3]

Collecting graphs by numerically integrating Eq. 11 of the paper;

ClearAll["Global`*"];
σ=1;
x0=1;
d=1;
k0=1;
h=-.4;
y=-1.7;
ys=0;
a=k0/(2d);
Table[1+a/2*Exp[a^2*(σ^2+d*t)]*NIntegrate[Exp[a*(Abs[x]+x0)]*(-Sign[(Abs[x]+x0)/(2(σ^2+d*t))+a]+Erf[(Abs[x]+x0+2*a*(σ^2+d*t))/(2*Sqrt[σ^2+d*t])]),{x,-500,500},Method->{"LocalAdaptive",Method->"LevinRule"},MaxRecursion->100,AccuracyGoal->10],{t,0,80}]
p1=ListPlot[%,PlotJoined->True,PlotLegends->Placed[{Style["(Subscript[k, 0])=1 nm/ps",13]},{Scaled[{0,0.5}],{y-0.2,h-2}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Dashed}];
Clear[k0];
k0=2;
p2=Plot[{1-a/2*Exp[a^2*(σ^2+d*t)]*NIntegrate[Exp[a*(Abs[x]+x0)]*Erfc[(Abs[x]+x0+2*a*(σ^2+d*t))/(2*Sqrt[σ^2+d*t])],{x,-1000,1000}]},{t,0,80},PlotLegends->Placed[{Style["(Subscript[k, 0])=2 nm/ps",13]},{Scaled[{0,0.5}],{y+0.1,h-1}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Thick}];
p3=Plot[{1-a/2*Exp[a^2*(σ^2+d*t)]*NIntegrate[Exp[a*(Abs[x]+x0)]*Erfc[(Abs[x]+x0+2*a*(σ^2+d*t))/(2*Sqrt[σ^2+d*t])],{x,-1000,1000}]},{t,0,80},PlotLegends->Placed[{Style["(Subscript[k, 0])=3 nm/ps",13]},{Scaled[{0,0.5}],{y+0.1,h}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Dotted}];
Show[p1,p2,p3]
Export["reactn.eps"]

Using the analytical equation Eq. 12 from the paper; We see some oscillations in the survival probability which is analyzed in the upcoming sections;

ClearAll["Global`*"];
σ=.5;
x0=3;
d=.5;
h=-.4;
y=-1.7;
ys=0;
a=k0/(2d);
P2[t_,k0_]=1+Exp[a*x0+a^2(σ^2+d*t)](Sign[x0/(2(σ^2+d*t))+a]-Erf[(x0+2a*(σ^2+d*t))/(2Sqrt[(σ^2+d*t)])])+2*HeavisideTheta[-x0-2*a(σ^2+d*t)]*Exp[-a^2(σ^2+d*t)]-Erfc[x0/(Sqrt[4(σ^2+d*t)])];
p1=Plot[{P2[t,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=1 nm/ps",13]},{Scaled[{0,0.5}],{y+2*0.1,h-1}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Dashed}];
p2=Plot[{P2[t,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=2 nm/ps",13]},{Scaled[{0,0.5}],{y+2*0.1,h}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Gray}];
p3=Plot[{P2[t,3]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=3 nm/ps",13]},{Scaled[{0,0.5}],{y+2*0.1,h+1}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Thick}];
Show[p1,p2,p3]
Export["reactndisc.eps",%]

Validating the analytical survival probability Eq. 12 using the numerical integration

Clear[σ,x0,d,a,k0];
σ=0.053;
x0=1;
d=1;
a=k0/(2d);
k0=1;
p1=Plot[{{1-a/2*Exp[a^2*(σ^2+d*t)]*NIntegrate[Exp[a*(Abs[x]+x0)]*Erfc[(Abs[x]+x0+2*a*(σ^2+d*t))/(2*Sqrt[σ^2+d*t])],{x,-1000,1000},MinRecursion->9]},{Exp[a*x]}},{t,0,100},PlotRange->{0,1},PlotTheme->{"Scientific","GrayColor"},PlotLabel->{"Survival probability vs time"},AxesLabel->{"time (ps)","Survival Probability P(t)"},AxesStyle->Thick,AxesStyle->Thick,PlotStyle->{Blue,Dashed},PlotLegends->{"Our solution"}]
Show[p1,p[t]]

Comparison of the solution with Green’s function method;
 sub task(a)- fit the initial Gaussian distribution using 3 dirac deltas
 subtask (b)- Use the 3 delta solutions in addition to compare with the Gaussian solution

ClearAll["Global`*"];
G[x_,y_,t_]=-k0/4*Exp[k0/2(Abs[x]+Abs[y])+k0^2*t/4]*Erfc[(k0*t+Abs[x]+Abs[y])/(2*Sqrt[t])]+1/(2*Sqrt[Pi*t])*Exp[-(x+y)^2/(t*4)]
α=100;
S[z_]=Sqrt[α/π]Exp[-α*z^2];
Integrate[S[z],{z,-Infinity,Infinity}]
Plot[S[z],{z,-6,6},PlotRange->Full];
NIntegrate[S[x]*Cos[x],{x,-10,10}]
ClearAll[k0,k1,k2,k3,σ];
ϵ=0.03;
b=-1;
a=-ϵ+b;
c=ϵ+b;
σ=0.053;
h=-0.2;
y=-1.5;
ys=0;
p1=Exp[-(z+1)^2/(4*σ^2)]/(4*Pi*σ^2)^(1/2);
Plot[p1,{z,-8,8},PlotRange->All];
Integrate[p1,{z,-Infinity,Infinity}]
k3=NIntegrate[p1,{z,-3,a}]
k1=NIntegrate[p1,{z,a,c}]
k2=NIntegrate[p1,{z,c,3}]
k0=1
P[x_,0]=k3*S[z-a]+k1*S[z-b]+k2*S[z-c];
Plot[{{P[x,0]},{p1}},{z,-1.4,-0.6},PlotRange->All]
Plot[{{k3*S[z-a]},{k1*S[z-b]},{k2*S[z-c]},{Exp[-(z-0)^2/(4*σ^2)]/(4*Pi*σ^2)}},{z,-1.4,-0.6},PlotRange->All];
G[x,a,t]
p2=Plot[NIntegrate[k3*G[x,a,t]+k1*G[x,a,t]+k2*G[x,c,t],{x,-30,30}],{t,0,100},PlotLegends->Placed[{Style["Numerical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Dashed}]
Clear[σ,x0,d,a,k0];
σ=0.053;
x0=1;
d=1;
a=k0/(2d);
k0=1;
p1=Plot[{1-a/2*Exp[a^2*(σ^2+d*t)]*NIntegrate[Exp[a*(Abs[x]+x0)]*Erfc[(Abs[x]+x0+2*a*(σ^2+d*t))/(2*Sqrt[σ^2+d*t])],{x,-1000,1000},MinRecursion->9]},{t,0,100},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{0,1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->Gray]
Show[p2,p1]
Export["numer1.eps",%,ImageResolution->1400,"CompressionLevel"->0]

THis section is for checking where from the oscillations in P[t] comes from;

ClearAll["Global`*"];
σ=0.5;
x0=1;
d=1;
k0=1;
a=k0/(2d);
P2[t_]=1+Exp[a*x0+a^2(σ^2+d*t)](Sign[x0/(2(σ^2+d*t))+a]-Erf[(x0+2a*(σ^2+d*t))/(2Sqrt[(σ^2+d*t)])])+2*HeavisideTheta[-x0-2*a(σ^2+d*t)]*Exp[-a^2(σ^2+d*t)]-Erfc[x0/(Sqrt[4(σ^2+d*t)])];
Plot[P2[t],{t,0,1000},PlotRange->Full]

Print["Secon term"]
Plot[Exp[a*x0+a^2(σ^2+d*t)](Sign[x0/(2(σ^2+d*t))+a]-Erf[(x0+2a*(σ^2+d*t))/(2Sqrt[(σ^2+d*t)])]),{t,0,1000},PlotRange->All]
Print["3 term"];Plot[2*HeavisideTheta[-x0-2*a(σ^2+d*t)]*Exp[-a^2(σ^2+d*t)],{t,0,100},PlotRange->All]
Print["4 term"];
Plot[Erfc[x0/(Sqrt[4(σ^2+d*t)])],{t,0,1000},PlotRange->All]

COmparison with the method of images solutions in limit ϵ0 -> ∞; Creating data for x0=10; similarly other graphs

ClearAll[B,A,x]
P[t_]=Limit[Erf[Abs[x]*Exp[-B*t]/Sqrt[2d/B(1-Exp[-2*B*t])]],B->0]
d=0.5
x0=10;
s0=1
b=10
x1=1
k=10000000
x0=10
Plot[{1-Erfc[x0/(Sqrt[4(s0^2+d*t)])]*Sum[(2*s-1)!!(-1)^s/(2^s(s0^2+d*t)*k^{2*s}),{s,0,Infinity}],P[t]},{t,0,10}]
t1=Column[Table[1-Erfc[x0/(Sqrt[4(s0^2+d*t)])],{t,0,50}]]
Export["x0=10.dat",t1]

Drawing the contour of the problem