Exact diffusion dynamics of a distribution in the piece-wise linear potential

Abstract:  We give time-domain calculation  of P(x,t) for the  Smoluchowski equation with piece-wise linear potential. We consider an absorbing  boundary at arbitrary  point  and finite strength. This problem has  been unsolved in the time-domain even for the delta function initial  condition.   The proposed new method overcomes the difficulties of taking inverse Laplace transform as well as  overcomes  the   difficulty of treating transcendental poles analytically. The exact analytic form of the time-dependent probability distribution is presented, and the results are compared with the Laplace domain solution.

Introduction:
Diffusion dynamics of a  probability distribution over a potential surface can mimic various processes related to chemistry , biology , physics etc. The moelcular level reactions are usually understood by the  motion of the reacting species along the reaction coordinate. In case of adiabatic picture,  the effect of metastability on the reaction can be probed whereas the effect of inter system crossings can be studied through the diabatic picture. So far in the literature, these effects are studied in the kinetic picture in deriving the rate constants of such reaction-diffusion systems.  There has been no exact  time-domain solution to investigate the dependence of parameters affecting the time-dependent survival probability. However, the problems are following reasons: i) The Smoluchowski equation governing the reaction-diffusion is usually solved using Laplace transform. And the inverse Laplace transform has been not feasible when one adds even a single boundary to the equation.  ii) For the models incorporating the metastability, the models are snot  solvable owing to the complicated form of the potentials. The processes are studied through the following one-dimensional equation written as,

((∂P (x, t))/(∂t))  =D(((∂^2)P)/(∂(x^2)))+(1/((Subscript[k, B])T)) (∂/(∂x))(P((∂U)/(∂x)))-(Subscript[k, 0])S(x)P(x,t),---(1)

where `D' is the diffusion coefficient, `x' is the reaction coordinate,  (Subscript[k, B]) is the Boltzmann constant, U(x) is the shape of the potential well, and depending the problem suitable sink function can be chosen. In this paper, we give exact time-domain method to explore the models of each of one type.  i)   model A, U(x)=b|x| with S(x)=δ(x-a), the particle moves along the piece-wise linear potential with a finite absorbing boundary placed at `a ’ . The model B is given by, ii) U(x)=b|x+a|θ(-x)+θ(x)b(a-x) with S(x)=0.  The model is construction  of piece-wise linear potentials with a metastable barrier separating two minima. as considered in [].The analytical solution for the above models involve complex inverse Laplace transform as well as difficulty in treating the transcendental poles. This paper presents a new time-domain method to calculate the exact time-dependent probability distribution for both  the models A and B, The proposed time-domain method uses random variable transformations to  flat potential  to derive the propagator. The solution is compared with the numerical results and also with the  Laplace domain solution [].
Let us consider the model A, which is given by the following equations,

((∂P (x, t))/(∂t))  =D(((∂^2)P)/(∂(x^2)))+(1/((Subscript[k, B])T)) (∂/(∂x))(P((∂|x|)/(∂x)))---(2 a)
															P(x,0)=δ(x+(Subscript[x, 0])),---(2 b)

Figure 1: The Schematic representation of the models A and B respectively.

Let us consider a simple diffusion equation,

((∂(OverTilde[P])((OverTilde[x]), (OverTilde[t])))/(Overscript[∂t,     ~]))  =D(((∂^2)(OverTilde[P]))/(∂((OverTilde[x])^2)))---(3 a)
															(OverTilde[P])((OverTilde[x]), 0)=δ((OverTilde[x])+(Subscript[x, 0]))(e^(g((Subscript[x, 0])sgn((Subscript[x, 0]))))),---(3 b)

for which the solution is known. After taking the transformation given by the following, The above equation is transformed into variables as given by the following,

(OverTilde[P])((OverTilde[x]),(OverTilde[t]))==v(x,t) (Subscript[P, 0])(x((OverTilde[x])),t((OverTilde[t])))((OverTilde[t])==(t/2)) ((OverTilde[x])==x) 
(v(x,t)==(e^(-(1/2) D (g^2) t-g(x sgn(x)))))----(4)

the Eq.3 becomes,

(1/2)DDerivative[2,0][(Subscript[P, 0])](x,t)+DgDerivative[1,0][(Subscript[P, 0])](x,t)==Derivative[0,1][(Subscript[P, 0])](x,t),x>0----5a
(1/2)DDerivative[2,0][(Subscript[P, 0])](x,t)-DgDerivative[1,0][(Subscript[P, 0])](x,t)==Derivative[0,1][(Subscript[P, 0])](x,t),x<0---5b

Here is the solubility test, initial condition test.

ClearAll["Global`*"];
σ=0;
FX[x1_,t_]=Assuming[x≠0,Integrate[(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2,{x,-∞,x1}]];
P[x1_,t_]=D[FX[x1,t],x1]
Print["Analytical check for the total solution in regions x>0, x<0 and x==0"]
Assuming[x≠0&&x>0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
Assuming[x≠0&&x<0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
Limit[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t],x->0]
Print["Initial check"]
x=-x0;
P[x,t]

Which is the equation for piece-wise linear potential centered at 0. The solution for the Eq. 3a for the initial condition given by Eq. 3b is given by,

(Subscript[P, 0])(x,t)==((e^(-(1/2) D (g^2) t-((((x+(Subscript[x, 0])))^2)/(2 D t))-g(x sgn(x)-(Subscript[x, 0])sgn((Subscript[x, 0])))))/((Sqrt[2 π]) (Sqrt[D t])))----(6)

As the random variable transformation from x to (OverTilde[x]) is linear, the probability distribution of the new variable can be obtained by mere substitution. It is seen that the obtained solution obeys Eq. 6 and the initial condition. Let us compare the obtained solution to the already available solution derived by VM Kenkre[] (Analysis of Confined Random, JPC, 2016);

ClearAll["Global`*"];
Γ= g d;
l=d/(2Γ);
σ=0;
u1[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0*RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))];
u2[x_,t_]=((Exp[-x RealSign[x]/l])/(4l))*Erfc[((x*RealSign[x]+x0*RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2) t]))];
u[x_,t_]=u1[x,t]+u2[x,t];
Print["Analytical check of VM kenkre solution"]
Assuming[x≠0,Simplify[d/2 D[u[x,t],{x,2}]+d g* D[u[x,t]*RealSign[x],x]-D[u[x,t],t]]]
P[x1_,t_]=(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Print["Analytical check of Eq. 6"]
Assuming[x≠ 0&&(x≥ 0&&x0≥0||x≤0&&x0≤0),Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
σ=0;
Print["The analytical difference is"]
Simplify[u1[x,t]-P[x,t],Assumptions->Element["Global`*",Reals]&&Element[g,Reals]&&g>0&&Element[d,Reals]&&d>0&&Element[t,Reals]&&t>0&&x>0&&x0>0]
g=1;
d=.1;
x0=5;
a=1;
Print["The analytical difference is"]
f[t_]=Simplify[u1[x,t]-P[x,t],Assumptions->Element["Global`*",Reals]&&Element[g,Reals]&&g>0&&Element[d,Reals]&&d>0&&Element[t,Reals]&&t>0&&x>0&&x0>0]
Plot[f[5],{x,-10,10},PlotRange->Full,PlotLegends->Automatic];
Plot[{u[x,100],P[x,0100]},{x,-10,10},PlotRange->Full,PlotLegends->Automatic];
p1=Plot[{u[x,.1],u[x,1],u[x,50]},{x,-10,10},PlotRange->Full,PlotLegends->Automatic];
p2=Plot[{P[x,.1],P[x,1],P[x,50],g Abs[x]},{x,-10,10},PlotRange->Full,PlotLegends->Automatic];
Print["graphical comparison between the solutions at different times"]
Show[p1,p2]

It is observed that the our solution is same except for a missing term. When comparing the solutions graphically, there are no differences in short times. For longer time, when the distribution reached the minimum of the potential, our solution differs from that of VM Kenkre. The process of “shape-shifting” is not fully incorporated. The shape-shifting occurs from the artifact of the following boundary condition, that is derivable from the equation.

(lim_(ϵ->0))(Subscript[(([D g sgn(x) P(x,t)+(1/2) D Derivative[1,0][P](x,t)])^(0+ϵ)), 0-ϵ])=0---(7)

Checking the solution given by Eq. 6 satisfies the above boundary condition or not in the following section

ClearAll["Global`*"];
Γ= g d;
l=d/(2Γ);
σ=0;
u[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))];
uplus[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))]/.
{x*RealSign[x]->x};
uminus[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))]/.{x*RealSign[x]->-x};
(*Print["Analytical check of kenkre"]
Assuming[x≠a,Simplify[d/2 D[u[x,t],{x,2}]+d g* D[u[x,t]*RealSign[x-a],x]-D[u[x,t],t]]]*)
Print["Continuity check"];
Simplify[Limit[uplus[x,t]-uminus[x,t],x->0]];
Print["Initial condition check"];
x=-x0;
u[x,t];
Clear[x];
Print["Boundary condition check of kenkre through 2 methods"]
Simplify[Limit[d/2 D[uplus[x,t],{x,1}]+d g*uplus[x,t]-d/2 D[uminus[x,t],{x,1}]-d g*uminus[x,t]*(-1),x->0]]
Simplify[d/2 D[uplus[x,t],{x,1}]+d g*uplus[x,t]-d/2 D[uminus[x,t],{x,1}]-d g*uminus[x,t]*(-1)/.x->0]
P[x_,t_]=(Exp[-g*((x*RealSign[x])-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Pplus[x_,t_]=(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.{t->t/2,(x*RealSign[x])->x};
Pminus[x_,t_]=(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.{t->t/2,(x*RealSign[x])->-x};
Print["Analytical check of Mr. R Saravanan"];
Assuming[x≠ 0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]];
σ=0;
Print["Continuity check"];
Simplify[Pplus[0,t]-Pminus[0,t]];
Print["Boundary check of Mr. R Saravanan which leads to differential equation at x=0"]
Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)-d g Pplus[0,t]/.x->0]
(*Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)/.x->a]
Simplify[-d g Pplus[a,t]/.x->a]
Simplify[-d g P[a,t]/.x->a]*)
Print["The analytical difference is"]
Simplify[u[x,t]-P[x,t],Assumptions->Element["Global`*",Reals]&&Element[g,Reals]&&g>0&&Element[d,Reals]&&d>0&&Element[t,Reals]&&t>0]
g=1;
d=.1;
x0=5;
a=-2;
Plot[{u[x,50],P[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic];
p1=Plot[{u[x,.1],u[x,1],u[x,10],u[x,500]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic];
p2=Plot[{P[x,.1],P[x,1],P[x,10],P[x,500]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic];
Show[p1,p2]

From inspection, we find that our solution rather obeys a boundary condition given by,

(lim_(ϵ->0))(Subscript[(([D g sgn(x) (Subscript[P, 0])(x,t)+(1/2) D Derivative[1,0][(Subscript[P, 0])](x,t)])^(0+ϵ)), 0-ϵ])=D g (Subscript[P, 0])(0,t)---(8)

The attempt to address the missing term can be either to solve the P(x,t) of the above equation. Another way is to assume a function (Subscript[P, 1])(x,t) such that the total solution is given by,P(x,t)=(Subscript[P, 0])(x,t)+(Subscript[P, 1])(x,t)  such that the function (Subscript[P, 1])(x,t) nullifies the extra term which can be expressed as,

(lim_(ϵ->0))(Subscript[(([D g sgn(x) (Subscript[P, 1])(x,t)+(1/2) D Derivative[1,0][(Subscript[P, 1])](x,t)])^(0+ϵ)), 0-ϵ])=-D g (Subscript[P, 0])(0,t)---(9)

However solving the Eq. 8 is much simpler. Let us simplify the Eq. (8) further to get,

(lim_(ϵ->0))(Subscript[(([(1/2) D Derivative[1,0][(Subscript[P, 0])](x,t)])^(0+ϵ)), 0-ϵ])+D g [(Subscript[P, 0])(0+ϵ,t)+(Subscript[P, 0])(0-ϵ,t)]=D g (Subscript[P, 0])(0,t)---(10)

using continuity the above equation simplifies to the following

(lim_(ϵ->0))(Subscript[(([(1/2) D Derivative[1,0][(Subscript[P, 0])](x,t)])^(0+ϵ)), 0-ϵ])-D g (Subscript[P, 0])(0,t)=0---(10)

The above solution can also be derived alternatively from the following equation,

d g (DiracDeltaSeq[x]) P(x,t)+(1/2) d Derivative[2,0][P](x,t)==Derivative[0,1][P](x,t)

For which the solution is readily available, adding this extra term to verify that the solution passes all the tests. The upcoming cells present the soluition of the above PDE in the transformed form

ClearAll["Global`*"];
σ=0;
eqn=(d/2)D[P[x,t],{x,2}]+d g P[x,t]DiracDelta[x]==D[P[x,t],t]//TraditionalForm(*eq.1*)
P[x,0]=Exp[g*RealAbs[x0]-g^2d t/2]DiracDelta[x-(Subscript[x, 0])]//TraditionalForm
eqn=(d/2)D[(Subscript[P, 1])[x,t],{x,2}]+2d g (Subscript[P, 1])[x,t]DiracDelta[x]+d g (Subscript[P, 0])[x,t]DiracDelta[x]==D[(Subscript[P, 1])[x,t],t]//TraditionalForm(*eq.2*)
(Subscript[P, 1])[x,0]=Exp[g*RealAbs[x0]-g^2d t/2]DiracDelta[x-(Subscript[x, 0])]//TraditionalForm
Print["For which we know the solution to be"](*multiplied Exp[g*(RealAbs[x0]-(x*RealSign[x-a]+2a*UnitStep[-x+a]))-g^2d t/2] to the soluton of eq. 1*)
ϵ0=((-d g)/(2(d/2)));
Pdif[x_,t_]=Simplify[((ϵ0)/2)Exp[g*(RealAbs[x0]-(x*RealSign[x]))-g^2d t/2]*Exp[ϵ0((x*RealSign[x])+RealAbs[x0])+ϵ0^2d/2t](-Sign[(1/(2d t/2))((x*RealSign[x])+RealAbs[x0])+ϵ0]+Erf[(x*RealSign[x]+RealAbs[x0]+2ϵ0 d/2 t)/(2Sqrt[d t/2])]),Assumptions->Element[d,Positive]&&Element[g,Positive]&&Element[t,Positive]&&Element[x0,Positive]]/.Sign[(1/(2d t/2))((x*RealSign[x])+RealAbs[x0])+ϵ0]->1
Print["An alternative solution using another PDE is"]
ϵ1=-2g;(*multiplied *Exp[2g RealAbs[x0]-g^2d t/2] to the soluton of eq. 2*)
Pdifanother[x_,t_]=Simplify[((- g)/2)Exp[2g RealAbs[x0]-4g^2d t/2]*Exp[ϵ1((x*RealSign[x])+RealAbs[x0])+ϵ1^2d/2t](-Sign[(1/(2d t/2))((x*RealSign[x])+RealAbs[x0])+ϵ1]+Erf[((x*RealSign[x])+RealAbs[x0]+2ϵ1 d/2 t)/(2Sqrt[d t/2])]),Assumptions->Element[d,Positive]&&Element[g,Positive]&&Element[t,Positive]&&Element[x0,Positive]]/.Sign[(1/(2d t/2))((x*RealSign[x])+RealAbs[x0])+ϵ1]->1
P0[x_,t_]=(Exp[-g*((x*RealSign[x])-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
P[x_,t_]=P0[x,t]+Pdif[x,t];
Pplus[x_,t_]=P[x,t]/.{x*RealSign[x]->x}
Pminus[x_,t_]=P[x,t]/.{x*RealSign[x]->-x}
Print["Boundary check of the new solution"]
Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)/.x->0,Assumptions->x0∈Positive]
Γ= g d;
l=d/(2Γ);
σ=0;
Print["The extra solution is compared with kenkre's one"]
udif[x_,t_]=(1/2) (E^(-2 g x RealSign[x])) g Erfc[((-d g t+x RealSign[x]+x0 RealSign[x0])/((Sqrt[2]) (Sqrt[d t])))]
(*((Exp[-RealAbs[x]/l])/(4l))*Erfc[((RealAbs[x]+RealAbs[x0]-Γ t)/(Sqrt[4 (σ^2+d/2) t]))]*)
g=1;
d=.1;
x0=-1;
a=-2;
Print["Checking the equality graphically"]
Plot[{udif[x,100],Pdif[x,100]},{x,-40,40},PlotRange->Full]
t={80};
Plot[Sign[-g+((RealAbs[x]+RealAbs[x0])/(d t))],{x,-40,40}]
(*Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0]+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[(1/(2Sqrt[(σ^2+d t/2)]))(x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0]+2b(σ^2+d t/2))]))/.RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1*)

Similarly the the piecewise linear potential with the sink at origin can be written by,

d g (DiracDeltaSeq[x]) P(x,t)-(Subscript[k, 0])P(x,t)+(1/2) d Derivative[2,0][P](x,t)==Derivative[0,1][P](x,t)

ClearAll["Global`*"];
σ=0;
Print["The solution is"]
b=(( k0-d g)/(2(d/2)));
P3[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(x*RealSign[x]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((x*RealSign[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((x*RealSign[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]))/.RealSign[((x*RealSign[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1;
P3[x,t]
Print["Solubility check of P3"]
Assuming[x≠ 0,Simplify[d/2 D[P3[x,t],{x,2}]+d g* D[P3[x,t]*RealSign[x],x]-D[P3[x,t],t]]]
P3plus[x_,t_]=P3[x,t]/.{x*RealSign[x]->x}
P3minus[x_,t_]=P3[x,t]/.{x*RealSign[x]->-x}
Print["Boundary check at x=0"]
Simplify[d/2 D[P3plus[x,t],{x,1}]+d g*P3plus[x,t]-d/2 D[P3minus[x,t],{x,1}]-d g*P3minus[x,t]*(-1)-k0*P3plus[x,t]/.x->0]
Print["The reduction test at (Subscript[k, 0])=0"]
k0=0;
Print["Boundary check with k0=0"]
Simplify[d/2 D[P3plus[x,t],{x,1}]+d g*P3plus[x,t]-d/2 D[P3minus[x,t],{x,1}]-d g*P3minus[x,t]*(-1)/.x->0]
a=0;
Simplify[Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])](0((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]))/.RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1]

Let us rederive the Swapan Gosh’s result of piecewise linear potential;

ClearAll["Global`*"]
y1[x_]=C1 Exp[(p+q)x];(*x<-b*)
(*y2[x_]=C3 Exp[(p+q)x]+C4 Exp[(q-p)x];(*0>x>-b*)*)
y3[x_]=A Exp[(p-q)x]+B*Exp[-(p+q)x];(*x0>x>0*)
{{C1,A,B}}={C1,A,B}/.Solve[{y1[0]==y3[0],
de(y3'[0]-y1'[0])+w de y3[0]-w de*y1[0](-1)==-1,de(-y3'[x0]-w y3[x0])==k0*y3[x0]},{C1,A,B}];
(*{{C1,C3,C4,A,B}}{C1,C3,C4,A,B}/.Solve[{y1[-b]==y2[-b],de(y2'[-b]-y1'[-b])==0,y2[0]==y3[0],de(y3'[0]-y2'[0])+w de y3[0]-w de*y2[0](-1)==-1,de(-y2'[x0]-w y2[x0])==k0*y2[x0]},{C1,C3,C4,A,B}]*)
L=(2p-w)(de(p-q)+(k0+w de))Exp[2 p x0]-w(de(p+q)-(k0+w de));
A2=((p+q-(k0+w de)/de)/L);
B2=((p-q+(k0+w de)/de)/L)Exp[2 p x0];
C2=A2+B2;
(*Simplify[C1];
Simplify[A];
Simplify[B];*)
q=w/2;
p=Sqrt[q^2+(s/de)];
Print["Check that the coefficients matches between both"]
Simplify[C1-C2]
Simplify[A-A2]
Simplify[B-B2]
y1a[x_]=C2 Exp[(p+q)x];(*x<0*)
y3a[x_]=A2 Exp[(p-q)x]+B2*Exp[-(p+q)x];(*x>0*)
Print["Boundary check of Mr Gosh"]
y1a[0]===y3a[0]
Simplify[de(y3a'[x0]+w y3a[x0])+k0*y3a[x0]]
Simplify[de(y3a'[0]-y1a'[0])+w de y3a[0]-w de*y1a[0](-1)+1]
ρ[s_]=(C2/(p+q))-(A2/(p-q))(1-Exp[(p-q)x0])+(B2/(p+q))(1-Exp[-(p+q)x0])
w=2;(*w=2g*)
de=.2/2;(*de=d/2*)
x0=0;(*position of the sink*)
k0=2;(*Sink strength*)
b=0;(*position of the initial condition*)
F[x_,s_]=Piecewise[{{y1[x],x<0},{y3[x],x0≥ x≥0}}];
f[x_,s_]=Integrate[F[x,s],x];
Q[s_]=FullSimplify[f[x0,s]];
(*Print["Check between the answers of survival prob"]
Plot[{ρ[s],Q[s],1/s},{s,1,20},PlotLegends->Automatic]*)
Attributes[NInverseLaplaceTransform]={NumericFunction};

NInverseLaplaceTransform:: wrongexpr="The argument `1` is not Numeric at fixed `2` .";NInverseLaplaceTransform:: wrongtime="The argument `1` is not Numeric.";
NInverseLaplaceTransform:: failed=
"The inverse could not be calculated. Suspect one of the following reasons: \n- the inverse does not exist; \n- the inverse or its derivatives have a discontinuity and the specified time falls near to it; \n- the inverse is zero; \n- the specified time is too large or small when the inverse is periodic-like.";

Options[NInverseLaplaceTransform]={"Startm"->5,"Method"->"Automatic"};

NInverseLaplaceTransform[expr_,s_Symbol,t_?NumericQ, opts:OptionsPattern[]]:=
Module[{F,z, x,M, f,g,fg,gg,ff,fgm, f12,g12, m,mmin,mmax,mb,meth,res }, 

mmin=Min[15,Max[5,OptionValue["Startm"]]];
mmax=Min[20,mmin+5];

meth=OptionValue["Method"];
If[Not[Or[meth==="FT",meth==="GWR"]],meth="Automatic"];

F[z_]=Rationalize[expr,0]/.s->z;
x=Rationalize[t,0];
If[Not[Block[{$MaxExtraPrecision=5000.},NumericQ[F[100]]]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongexpr,expr,s]);res=$Failed;Return[res]];If[Not[And[NumericQ[t],t>0]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongtime,t]);res=$Failed;Return[res]];
f12=g12=Array[0&,mmax];
Do[
M=2^m;mb=m;ff=fgm=fg=gg=False;

If[Or[meth==="FT",meth==="Automatic"],
f=FT[F, x, M]; 
f12[[m]]=If[$MinMachineNumber<Abs[f]<$MaxMachineNumber,f//N,0.];
ff=And[f12[[m]]≠0,f12[[m]]==f12[[m-1]]];
fgm=And[f12[[m]]≠0,f12[[m]]==g12[[m-1]]];
     If[Or[ff,fgm],Break[]];      
     If[And[m==mmax-1,meth==="FT",f12[[m]]==0],res=$Failed;Break[]];
     If[And[m==mmax,f12[[m]]==0,g12[[m-1]]==0],res=$Failed;Break[]];
];

If[Or[meth==="GWR",meth==="Automatic"],
g=GWR[F, x, M];
g12[[m]]=If[$MinMachineNumber<Abs[g]<$MaxMachineNumber,g//N,0.];
fg=And[g12[[m]]≠0,f12[[m]]==g12[[m]]];
gg=And[g12[[m]]≠0,g12[[m]]==g12[[m-1]]];
      If[Or[fg,gg],Break[]];
      If[And[m==mmax-1,meth==="GWR",g12[[m]]==0],res=$Failed;Break[]];   
      If[And[m==mmax,g12[[m]]==0],res=$Failed;Break[]];
];
 If[m==mmax,res=$Failed]
,{m,mmin,mmax}
];
If[Or[fg,gg],res=g12[[mb]]; Return[res]];
If[Or[ff,fgm],res=f12[[mb]]; Return[res]];
res==$Failed;
(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform::failed]);
res
];

GWR[F_, t_?NumericQ, M_:32, precin_:0] := Module[
{ fv,w, wm,wp, best, tau = Log[2]/t,  prec,k,n,i}, 
If[precin ≤ 0, prec = 2.1*M, prec = precin]; 
If[prec ≤ $MachinePrecision, prec = $MachinePrecision]; 
If[Precision[tau] < prec, tau = SetPrecision[tau, prec]]; 
Block[{$MaxExtraPrecision=5000.},
fv=Table[ N[F[i tau], prec], {i, 2M}]; 
w=Table[ 
tau*((2*n)!/(n!*(n - 1)!))*Total[Table[Binomial[n, i]*(-1)^i*fv[[n + i]], {i, 0, n}]]
, {n,  M}]; 
];
best=Last[w];
If[Not[NumberQ[best]],best=$Failed;Return[best]];
wm=Table[0,{M}];
Do[
wp=Quiet[k/Differences[w]];
wp=Quiet[Drop[wm,1]+wp];
wm=Drop[w,-1];
w=wp;
If[EvenQ[k],If[NumberQ[Last[w]],best =Last[w],Break[],Break[]]]
,{k,  M -1}]; 
If[Not[NumberQ[best]],best=$Failed];
best
]

FT[F_,t_?NumericQ,M_:32]:=Module[{np,p,r,s,theta,sigma,tab,f},
np=Max[M,$MachinePrecision];
p=SetPrecision[t,np];
Block[{$MaxExtraPrecision=5000.},
r=2M/(5p);
s=r theta (Cot[theta]+I);
sigma=theta+(theta Cot[theta]-1)Cot[theta];
Quiet[
tab=Flatten[{Exp[p r] F[r]/2,Table[Exp[p s](1+I sigma)F[s],{theta,Pi/M,(M-1)Pi/M,Pi/M}]/.ComplexInfinity->0}];
f=Total[Re[(r/M) tab]]
]
];
f
]
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}]
ab2=ListPlot[%,PlotLegends->"Goshmy",PlotStyle->Green]
Table[NInverseLaplaceTransform[ρ[s],s,t],{t,.1,50}]
ab1=ListPlot[%,PlotLegends->"Gosh him",PlotStyle->Red]

(Legended[, "Gosh him"])

NOw let us the add the effect of arbitrarily placed initial condition; (NOTE: It matches the above solution under limits and all boundary cnditions and initial test)

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
y1[x_]=C1 Exp[(p+q)x];(*x<-b*)
y2[x_]=C3 Exp[(p+q)x]+C4 Exp[(q-p)x];(*0>x>-b*)
y3[x_]=A Exp[(p-q)x]+B*Exp[-(p+q)x];(*x0>x>0*)
(*{{C1,A,B}}={C1,A,B}/.Solve[{y1[0]==y3[0],
de(y3'[0]-y1'[0])+w de y3[0]-w de*y1[0](-1)==-1,de(-y3'[x0]-w y3[x0])==k0*y3[x0]},{C1,A,B}];*)
{{C1,C3,C4,A,B}}={C1,C3,C4,A,B}/.Solve[{y1[-b]==y2[-b],de(y2'[-b]-y1'[-b])==-1,y2[0]==y3[0],de(y3'[0]-y2'[0])+w de y3[0]-w de*y2[0](-1)==0,de(-y3'[x0]-w y3[x0])==k0*y3[x0]},{C1,C3,C4,A,B}];
q=w/2;
p=Sqrt[q^2+(s/de)];
w=2g;(*w=2g*)
de=.2/2;(*de=d/2*)
x0=0;(*position of the sink*)
k0=3;(*Sink strength*)
Clear[b];
b=5;(*position of the initial condition*)
F[x_,s_]=Piecewise[{{y1[x],x<-b},{y2[x],0≥ x≥-b},{y3[x],x0≥ x≥0}}];
f[x_,s_]=Integrate[F[x,s],x];
Q[s_]=FullSimplify[f[x0,s]];
(*Print["Check between the answers of survival prob"]
Plot[{ρ[s],Q[s],1/s},{s,1,20},PlotLegends->Automatic]*)
Attributes[NInverseLaplaceTransform]={NumericFunction};

NInverseLaplaceTransform:: wrongexpr="The argument `1` is not Numeric at fixed `2` .";NInverseLaplaceTransform:: wrongtime="The argument `1` is not Numeric.";
NInverseLaplaceTransform:: failed=
"The inverse could not be calculated. Suspect one of the following reasons: \n- the inverse does not exist; \n- the inverse or its derivatives have a discontinuity and the specified time falls near to it; \n- the inverse is zero; \n- the specified time is too large or small when the inverse is periodic-like.";

Options[NInverseLaplaceTransform]={"Startm"->5,"Method"->"Automatic"};

NInverseLaplaceTransform[expr_,s_Symbol,t_?NumericQ, opts:OptionsPattern[]]:=
Module[{F,z, x,M, f,g,fg,gg,ff,fgm, f12,g12, m,mmin,mmax,mb,meth,res }, 

mmin=Min[15,Max[5,OptionValue["Startm"]]];
mmax=Min[20,mmin+5];

meth=OptionValue["Method"];
If[Not[Or[meth==="FT",meth==="GWR"]],meth="Automatic"];

F[z_]=Rationalize[expr,0]/.s->z;
x=Rationalize[t,0];
If[Not[Block[{$MaxExtraPrecision=5000.},NumericQ[F[100]]]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongexpr,expr,s]);res=$Failed;Return[res]];If[Not[And[NumericQ[t],t>0]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongtime,t]);res=$Failed;Return[res]];
f12=g12=Array[0&,mmax];
Do[
M=2^m;mb=m;ff=fgm=fg=gg=False;

If[Or[meth==="FT",meth==="Automatic"],
f=FT[F, x, M]; 
f12[[m]]=If[$MinMachineNumber<Abs[f]<$MaxMachineNumber,f//N,0.];
ff=And[f12[[m]]≠0,f12[[m]]==f12[[m-1]]];
fgm=And[f12[[m]]≠0,f12[[m]]==g12[[m-1]]];
     If[Or[ff,fgm],Break[]];      
     If[And[m==mmax-1,meth==="FT",f12[[m]]==0],res=$Failed;Break[]];
     If[And[m==mmax,f12[[m]]==0,g12[[m-1]]==0],res=$Failed;Break[]];
];

If[Or[meth==="GWR",meth==="Automatic"],
g=GWR[F, x, M];
g12[[m]]=If[$MinMachineNumber<Abs[g]<$MaxMachineNumber,g//N,0.];
fg=And[g12[[m]]≠0,f12[[m]]==g12[[m]]];
gg=And[g12[[m]]≠0,g12[[m]]==g12[[m-1]]];
      If[Or[fg,gg],Break[]];
      If[And[m==mmax-1,meth==="GWR",g12[[m]]==0],res=$Failed;Break[]];   
      If[And[m==mmax,g12[[m]]==0],res=$Failed;Break[]];
];
 If[m==mmax,res=$Failed]
,{m,mmin,mmax}
];
If[Or[fg,gg],res=g12[[mb]]; Return[res]];
If[Or[ff,fgm],res=f12[[mb]]; Return[res]];
res==$Failed;
(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform::failed]);
res
];

GWR[F_, t_?NumericQ, M_:32, precin_:0] := Module[
{ fv,w, wm,wp, best, tau = Log[2]/t,  prec,k,n,i}, 
If[precin ≤ 0, prec = 2.1*M, prec = precin]; 
If[prec ≤ $MachinePrecision, prec = $MachinePrecision]; 
If[Precision[tau] < prec, tau = SetPrecision[tau, prec]]; 
Block[{$MaxExtraPrecision=5000.},
fv=Table[ N[F[i tau], prec], {i, 2M}]; 
w=Table[ 
tau*((2*n)!/(n!*(n - 1)!))*Total[Table[Binomial[n, i]*(-1)^i*fv[[n + i]], {i, 0, n}]]
, {n,  M}]; 
];
best=Last[w];
If[Not[NumberQ[best]],best=$Failed;Return[best]];
wm=Table[0,{M}];
Do[
wp=Quiet[k/Differences[w]];
wp=Quiet[Drop[wm,1]+wp];
wm=Drop[w,-1];
w=wp;
If[EvenQ[k],If[NumberQ[Last[w]],best =Last[w],Break[],Break[]]]
,{k,  M -1}]; 
If[Not[NumberQ[best]],best=$Failed];
best
]

FT[F_,t_?NumericQ,M_:32]:=Module[{np,p,r,s,theta,sigma,tab,f},
np=Max[M,$MachinePrecision];
p=SetPrecision[t,np];
Block[{$MaxExtraPrecision=5000.},
r=2M/(5p);
s=r theta (Cot[theta]+I);
sigma=theta+(theta Cot[theta]-1)Cot[theta];
Quiet[
tab=Flatten[{Exp[p r] F[r]/2,Table[Exp[p s](1+I sigma)F[s],{theta,Pi/M,(M-1)Pi/M,Pi/M}]/.ComplexInfinity->0}];
f=Total[Re[(r/M) tab]]
]
];
f
]
h=-.4;
y=-1.1;
ys=0;
g=1
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}];
ab2=ListPlot[%,PlotStyle->Black,PlotLegends->Placed[{Style["g=0.5",13]},{Scaled[{0,0.5}],{y+1,h-1}}]];

We saw that the solution constructed for sink placed at the minima of the |x| potential, it satisfies all the solubility, initial condition test, and spatial boundary conditions. Let us compare with the numerical results for now. And then proceed to calculate for the arbitrarily placed sink. Following are the numerical tests: 1) Compare with Swapan Gosh inverse Laplace transform, 2) Compare with NDSolve.

Protect[ab2]
ClearAll["Global`*"]
σ=0;
b=(( k0-d g)/(2(d/2)));
P1[x_,t_]=Exp[-g*(-RealAbs[x0]+RealAbs[x])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))-(b/2)*Exp[b(RealAbs[x]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(Erfc[((RealAbs[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]));
P2[x_,t_]=Exp[-g*(-RealAbs[x0]+RealAbs[x])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(RealAbs[x]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((RealAbs[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((RealAbs[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]));
Qt=(1/(2 (Sqrt[d t])))(E^(g (-x0+Abs[x0]))) (((1+(E^(2 g x0)))/(Sqrt[(1/(d t))]))-(((d g t-x0) Erf[((Sqrt[((((-d g t+x0))^2)/(d t))])/(Sqrt[2]))])/(Sqrt[((((-d g t+x0))^2)/(d t))]))-(((E^(2 g x0)) (d g t+x0) Erf[((Sqrt[((((d g t+x0))^2)/(d t))])/(Sqrt[2]))])/(Sqrt[((((d g t+x0))^2)/(d t))])))+ϵ1*Exp[(ϵ1+g)RealAbs[x0]-d g^2/2t+ϵ1^2 d t/2]((1/(ϵ1-g)))(Sign[(x0/(d t))+ϵ1]-Erf[((x0+ϵ1 d t)/(Sqrt[2 d t]))]+2Exp[(ϵ1-g)(-x0-ϵ1 d t)]UnitStep[-x0-ϵ1 d t]-Exp[((g-ϵ1)/2)(2 x0+d t(g+ϵ1))]*Erfc[((x0+d g t)/(Sqrt[2d t]))]);
P[x,t];
g=1;
d=.2;
x0=5;
k0=3;
h=-.4;
y=-1.1;
ys=0;
ϵ1=b;
p2=Plot[Qt,{t,1,50},PlotRange->Full,PlotLegends->Placed[{Style["Q(t)",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}]];
p0=Plot[{P[x,.1],P[x,1],P[x,10],P[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic];
Flatten[Table[NIntegrate[P1[x,t],{x,-20,20},AccuracyGoal->Infinity,Method->"Oscillatory"],{t,1,50}]];
p1=ListPlot[%,PlotJoined->True,PlotRange->Full,PlotLegends->Placed[{Style["Qnc[t]",13]},{Scaled[{0,0.5}],{y-4.5ys,h+2}}]];
Flatten[Table[NIntegrate[P2[x,t],{x,-20,20},AccuracyGoal->Infinity,Method->"Oscillatory"],{t,1,50}]]
p4=ListPlot[%,PlotJoined->True,PlotRange->Full,PlotLegends->Placed[{Style["Q[t]NInt",13]},{Scaled[{0,0.5}],{y-4.5ys,h+3}}]]
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
dThick=40;
dTime=50;
eqSimplee=D[Pe[x,t],t]==d/2*D[Pe[x,t],{x,2}]+RealSign[x]*d g*D[Pe[x,t],x]+(2d g-k0)*S[x]*Pe[x,t];
eqInie=Pe[x,0]==√(1000/π)*Exp[-1000*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
h=-.4;
y=-1.5;
ys=0;
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqInie,eqBlefte,eqBleftg},{Pe[x,t]},{x,-dThick,dThick},{t,0,dTime}]
Plot3D[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick},{t,0,dTime},PlotRange->Full]
p3=Plot[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick}],{t,0,50},PlotLegends->Placed[{Style["NDSolve",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
Show[p3,p4,p1,p2,ab2]
Show[p3,p1,p2,ab2]
Beep[]

Flatten[Table[NIntegrate[P1[x,t],{x,-20,20},AccuracyGoal->Infinity,Method->"Oscillatory"],{t,1,50}]];
p1=ListPlot[%,PlotRange->Full,PlotLegends->Placed[{Style["Qnc[t]",13]},{Scaled[{0,0.5}],{y-4.5ys,h+2}}]];
Show[p3,p1,p2,ab2]

Show[p2,ab2,p1]

ListPlot[]

Show[ab1,p1]

Let us generate some useful graphs for the analysis

ClearAll["Global`*"]
σ=0;
b=(( k0-d g)/(2(d/2)));

ϵ1=b;
Qt[t_]=(1/(2 (Sqrt[d t])))(E^(g (-x0+Abs[x0]))) (((1+(E^(2 g x0)))/(Sqrt[(1/(d t))]))-(((d g t-x0) Erf[((Sqrt[((((-d g t+x0))^2)/(d t))])/(Sqrt[2]))])/(Sqrt[((((-d g t+x0))^2)/(d t))]))-(((E^(2 g x0)) (d g t+x0) Erf[((Sqrt[((((d g t+x0))^2)/(d t))])/(Sqrt[2]))])/(Sqrt[((((d g t+x0))^2)/(d t))])))+ϵ1*Exp[(ϵ1+g)RealAbs[x0]-d g^2/2t+ϵ1^2 d t/2]((1/(ϵ1-g)))(Sign[(x0/(d t))+ϵ1]-Erf[((x0+ϵ1 d t)/(Sqrt[2 d t]))]+2Exp[(ϵ1-g)(-x0-ϵ1 d t)]UnitStep[-x0-ϵ1 d t]-Exp[((g-ϵ1)/2)(2 x0+d t(g+ϵ1))]*Erfc[((x0+d g t)/(Sqrt[2d t]))]);
g=.1;
d=.5;
x0=5;
k0=3;
h=-.4;
y=-1.1;
ys=0;
p1=Plot[Qt[t],{t,0,50},PlotRange->Full,PlotStyle->Gray,PlotLegends->Placed[{Style["g=0.1",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotTheme->{"Scientific"}];
ClearAll[f,g]
g=.5;
p2=Plot[Qt[t],{t,0,50},PlotRange->Full,PlotStyle->Dashed,PlotLegends->Placed[{Style["g=0.5",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotTheme->"Scientific"];
ClearAll[f,g]
g=.8;
p3=Plot[Qt[t],{t,0,50},PlotRange->Full,PlotStyle->{Thick},PlotLegends->Placed[{Style["g=0.8",13]},{Scaled[{0,0.5}],{y-4.5ys,h}}],PlotTheme->{"Scientific","Monochrome"}];
ClearAll[f,g]
g=1;
p4=Plot[Qt[t],{t,0,50},PlotRange->Full,PlotStyle->{Dotted,Black},PlotLegends->Placed[{Style["g=1",13]},{Scaled[{0,0.5}],{y-4.5ys,h+1}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}];
Show[p4,p1,p2,p3]

There are some term missing from our derived solution. However that can be catched by using Laplace domain and comparing with ours. Now let us solve the difficult problem of arbitrarily placed sink.

d g δ(x)P(x,t)+d g (DiracDeltaSeq[x-a]) P(x,t)-(Subscript[k, 0])(DiracDeltaSeq[x-a])P(x,t)+(1/2) d Derivative[2,0][P](x,t)==Derivative[0,1][P](x,t)

THe model is similar to that of the Oster-Nishijima’s model. And the solution is already derived by us.

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
Print["Our solution is"]
σ=0;
ϵ1=(k1/(2(d/2)));
ϵ2=(k2/(2(d/2)));
k1=-d g;
k2=k0;
Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
Heaviside[x_]:=Piecewise[{{0,x≤ 0},{1,x>0}}];
P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](-RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(((x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t)))/(2Sqrt[(σ^2+ d t)]))])+Heaviside[a-x](((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((x+2(b-a)+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((-x+2a+x0)/(2Sqrt[(σ^2+ d t)]))]))/.{RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]->1,t->t/2};Pe[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x])-g^2d t/2]*P1[x,t]
Pebplus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]->x,Heaviside[a-x]->0,x*RealSign[x]->x};
Pebminus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->0,x*RealSign[x]->x};
Peaplus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->0,x*RealSign[x]->x};
Peaminus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->1,x*RealSign[x]->-x};
a=0;
b=xc;
Print["Solubility check of Pe"]
Assuming[x≠ a&&x≠b&&x0>0&&xc>0,Simplify[d/2 D[Pe[x,t],{x,2}]+d g* D[Pe[x,t]*RealSign[x],x]-D[Pe[x,t],t]]]
Print["Continuity check of P3 a, b"]
Simplify[Pebplus[b,t]-Pebminus[b,t]]
Simplify[Peaplus[a,t]-Peaminus[a,t]]
a=0;
b=xc;
Print["Boundary check for Pe"]
Assuming[b≥a,Simplify[d/2 D[Peaplus[x,t],{x,1}]+d g*Peaplus[x,t]-d/2 D[Peaminus[x,t],{x,1}]-d g*Peaplus[x,t]*(-1)/.x->a]]
Assuming[b≥a,Simplify[d/2 D[Pebplus[x,t],{x,1}]+d g*Pebplus[x,t]-d/2 D[Pebminus[x,t],{x,1}]-d g*Pebminus[x,t]-(k0)Pebplus[b,t]/.x->b]]
Beep[]

ClearAll["Global`*"]
σ=0;
P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))E^(-(((x+x0)^2)/(4 (σ^2+ d t))))+((ϵ2)/2)E^(ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0))(-RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(((x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t)))/(2Sqrt[(σ^2+ d t)]))])+Heaviside[a-x](((ϵ1)/2)E^(ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0))(-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)E^(ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0))(-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((x+2(b-a)+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((-x+2a+x0)/(2Sqrt[(σ^2+ d t)]))]))/.{t->t/2,d->D,ϵ2->(Subscript[ϵ, 2]),ϵ1->(Subscript[ϵ, 1]),Unit->U,x0->(Subscript[x, 0]),a->0,b->a,Heaviside->H};ϵ1=(k1/(2(D/2)));
ϵ2=(k2/(2(D/2)));
k1=-D g;
k2=(Subscript[k, 0]);
Pe[x_,t_]=E^(-g*(-RealAbs[x0]+x*RealSign[x])-g^2d t/2)*P1[x,t]
TeXForm[%]

We have got to check all the tests, Now numerical test

Protect[ab2];
ClearAll["Global`*"]
σ=0;
ϵ1=(k1/(2(d/2)));
ϵ2=(k2/(2(d/2)));
k1=-d g;
k2=k0;
Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
Heaviside[x_]:=Piecewise[{{0,x≤ 0},{1,x>0}}];
P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](-RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+Heaviside[a-x](((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((x+2(b-a)+x0))/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((-x+2a+x0))/(2Sqrt[(σ^2+ d t)]))]))/.{t->t/2};(*{RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]->1,t->t/2}*)Pnc[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]-((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](Erfc[(x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+Heaviside[a-x](((-ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](Erfc[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](Erfc[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((x+2(b-a)+x0))/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((-x+2a+x0))/(2Sqrt[(σ^2+ d t)]))]))/.{t->t/2};P[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x])-g^2d t/2]*Pnc[x,t];
Qt[t_]=((Exp[g(RealAbs[x0]-x0)])/2)(Erfc[((Sqrt[(-x0+d g t)^2])/(Sqrt[2 d t]))]+Exp[2 g x0]Erfc[((d g t+x0)/(Sqrt[2 d t]))])+((Exp[ϵ2^2d t/2+ϵ2 x0+g RealAbs[x0]-g^2 d t/2])/2)(ϵ2 Exp[ϵ2*2b](-RealSign[((2b+x0)/(d t))+ϵ2]+Erf[((2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])((1/(g-ϵ2))+(1/(g+ϵ2)))+ϵ2 Exp[(ϵ2-g)b](-RealSign[((b+x0)/(d t))+ϵ2]+Erf[((b+x0+ϵ2 d t)/(Sqrt[2 d t]))])(-(1/(ϵ2-g))-(1/(ϵ2+g)))-ϵ2 Exp[ϵ2(-x0-ϵ2 d t)]((2/(g-ϵ2))Exp[g(2b+x0+ϵ2 d t)]UnitStep[-2b-x0-ϵ2 d t]+(2/(-(ϵ2+g)))Exp[-g(2b+x0+ϵ2 d t)]UnitStep[2b+x0+ϵ2 d t]UnitStep[b-2b-x0-ϵ2 d t]-((2 Exp[g(x0+ϵ2 d t)])/(ϵ2-g))UnitStep[-x0-ϵ2 d t-b])+ϵ2(((1/(g-ϵ2)))Exp[2b g+((g-ϵ2)/2)(2x0+(g+ϵ2)d t)]*Erfc[((2b+x0+d g t)/(Sqrt[2 d t]))]-((Exp[((g-ϵ2)/2)(2x0+(g+ϵ2)d t)])/(ϵ2-g))*Erfc[((b+d g t+x0)/(Sqrt[2 d t]))]+((Exp[-2b g+(1/2)(-2x0+d t(g-ϵ2))(g+ϵ2)])/(g+ϵ2))(Erfc[((Sqrt[(2b-d g t+x0)^2])/(Sqrt[2 d t]))]-Erfc[((Sqrt[(b-d g t+x0)^2])/(Sqrt[2 d t]))]))+ϵ1*Exp[2b ϵ2](((Exp[-ϵ2 a+g a])/(g-ϵ2))(-Sign[((-a+2b+x0)/(d t))+ϵ2]+Erf[((-a+2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])-(2/(g-ϵ2))Exp[(g-ϵ2)(2b+x0+ϵ2 d t)]UnitStep[a-2b-x0-ϵ2 d t]+(2/(2(g-ϵ2)))Exp[((g-ϵ2)/2)(4b+2x0+d t(g+ϵ2))]*Erfc[((2b-a+d g t+x0)/(Sqrt[2d t]))])-ϵ1*Exp[ϵ2 2(b-a)](((Exp[(ϵ2+g) a])/(g+ϵ2))(-Sign[((a+2(b-a)+x0)/(d t))+ϵ2]+Erf[((a+2(b-a)+x0+ϵ2 d t)/(Sqrt[2 d t]))])-((Exp[-(g+ϵ2)(2(b-a)+x0+ϵ2 d t)])/(g+ϵ2))(-2)UnitStep[-a+2b+x0+ϵ2 d t]-((Exp[((g+ϵ2)/2)(4a-4b+d g t-x0-ϵ2 d t)])/(g+ϵ2))Erfc[(((a-2b+d g t-x0))/(Sqrt[2d t]))]));
Qtnc[t_]=((Exp[g(RealAbs[x0]-x0)])/2)(Erfc[((Sqrt[(-x0+d g t)^2])/(Sqrt[2 d t]))]+Exp[2 g x0]Erfc[((d g t+x0)/(Sqrt[2 d t]))])+((Exp[ϵ2^2d t/2+ϵ2 x0+g RealAbs[x0]-g^2 d t/2])/2)(ϵ2 Exp[ϵ2*2b](-Erfc[((2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])((1/(g-ϵ2))+(1/(g+ϵ2)))+ϵ2 Exp[(ϵ2-g)b](-Erfc[((b+x0+ϵ2 d t)/(Sqrt[2 d t]))])(-(1/(ϵ2-g))-(1/(ϵ2+g)))-ϵ2 Exp[ϵ2(-x0-ϵ2 d t)]((2/(g-ϵ2))Exp[g(2b+x0+ϵ2 d t)]UnitStep[-2b-x0-ϵ2 d t]+(2/(-(ϵ2+g)))Exp[-g(2b+x0+ϵ2 d t)]UnitStep[2b+x0+ϵ2 d t]UnitStep[b-2b-x0-ϵ2 d t]-((2 Exp[g(x0+ϵ2 d t)])/(ϵ2-g))UnitStep[-x0-ϵ2 d t-b])+ϵ2(((1/(g-ϵ2)))Exp[2b g+((g-ϵ2)/2)(2x0+(g+ϵ2)d t)]*Erfc[((2b+x0+d g t)/(Sqrt[2 d t]))]-((Exp[((g-ϵ2)/2)(2x0+(g+ϵ2)d t)])/(ϵ2-g))*Erfc[((b+d g t+x0)/(Sqrt[2 d t]))]+((Exp[-2b g+(1/2)(-2x0+d t(g-ϵ2))(g+ϵ2)])/(g+ϵ2))(Erfc[((Sqrt[(2b-d g t+x0)^2])/(Sqrt[2 d t]))]-Erfc[((Sqrt[(b-d g t+x0)^2])/(Sqrt[2 d t]))]))+0ϵ1*Exp[2b ϵ2](((Exp[-ϵ2 a+g a])/(g-ϵ2))(-Erfc[((-a+2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])-(2/(g-ϵ2))Exp[(g-ϵ2)(2b+x0+ϵ2 d t)]UnitStep[a-2b-x0-ϵ2 d t]+(2/(2(g-ϵ2)))Exp[((g-ϵ2)/2)(4b+2x0+d t(g+ϵ2))]*Erfc[((2b-a+d g t+x0)/(Sqrt[2d t]))])-ϵ1*Exp[ϵ2 2(b-a)](((Exp[(ϵ2+g) a])/(g+ϵ2))(-Erfc[((a+2(b-a)+x0+ϵ2 d t)/(Sqrt[2 d t]))])-((Exp[-(g+ϵ2)(2(b-a)+x0+ϵ2 d t)])/(g+ϵ2))(-2)UnitStep[-a+2b+x0+ϵ2 d t]-((Exp[((g+ϵ2)/2)(4a-4b+d g t-x0-ϵ2 d t)])/(g+ϵ2))Erfc[(((a-2b+d g t-x0))/(Sqrt[2d t]))]));
g=1;
d=.3;
x0=4;
k0=3;
a=0;
b=xc;
xc=-1;
h=-.4;
y=-1.1;
ys=0;
(*-((ϵ1)/2)(-RealSign[((x+x0)/(2(σ^2+d t)))]+Erf[((x+x0)/(2Sqrt[(σ^2+ d t)]))]-RealSign[((-x+2b+x0)/(2(σ^2+d t)))]+Erf[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]+RealSign[((x+2(b-a)+x0)/(2(σ^2+d t)))]-Erf[(((x+2(b-a)+x0))/(2Sqrt[(σ^2+ d t)]))]+RealSign[((-x+2a+x0)/(2(σ^2+d t)))]-Erfc[(((-x+2a+x0))/(2Sqrt[(σ^2+ d t)]))])*)
Ex[t_]:=ϵ1/2*NIntegrate[Erfc[((Sqrt[(x+x0)^2])/(Sqrt[2 d t]))]+Erfc[((Sqrt[(-x+2b+x0)^2])/(Sqrt[2d t]))]-Erfc[((Sqrt[(x+2(b-a)+x0)^2])/(Sqrt[2d t]))]-Erfc[((Sqrt[(-x+2a+x0)^2])/(Sqrt[2d t]))],{x,-∞,a},Method->"AdaptiveMonteCarlo",MaxRecursion->20]
c=f/.Flatten[Block[{t=0.1},Solve[f*(Qt[t]+Ex[t])==1,f]]];
Table[c(Qt[t]+Ex[t]),{t,0.1,50}]
p1=ListPlot[%,PlotJoined->True,PlotRange->Full,PlotStyle->{Red,Dashed},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotTheme->"Scientific",PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]}]
e=f/.Flatten[Block[{t=0.1},Solve[f*(Qtnc[t]+Ex[t])==1,f]]];
Table[e(Qtnc[t]+0Ex[t]),{t,0.1,50}]
p2=ListPlot[%,PlotRange->Full,PlotLegends->"nocorrection",PlotStyle->Brown]
p3=Plot[{P[x,.1],P[x,1],P[x,10],P[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic]
m=l/.Solve[l*NIntegrate[P[x,.15],{x,-30,30},AccuracyGoal->Infinity,Method->"AdaptiveQuasiMonteCarlo",MaxRecursion->20]==1,l](*"AdaptiveMonteCarlo"*)
Flatten[Table[m NIntegrate[ P[x,t],{x,-30,30},AccuracyGoal->Infinity,Method->"AdaptiveQuasiMonteCarlo",MaxRecursion->20],{t,.1,50}]]
p4=ListPlot[%,PlotJoined->True,PlotRange->Full,PlotLegends->Placed[{Style["NIntegrate",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}]]
Show[p4,ab2,p2,p1]
Show[p1,ab2]
(*S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
dThick=20;
dTime=50;
eqSimplee=D[Pe[x,t],t]==d/2*D[Pe[x,t],{x,2}]+d g RealSign[x]*D[Pe[x,t],x]+2d g S[x]*Pe[x,t]-k0*S[x-xc]*Pe[x,t];
eqInie=Pe[x,0]==√(1000/π)*Exp[-1000*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
h=-.4;
y=-1.5;
ys=0;
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqInie,eqBlefte,eqBleftg},{Pe[x,t]},{x,-dThick,dThick},{t,0,dTime}]
Table[{t,Flatten[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick},AccuracyGoal->10,PrecisionGoal->6,Method->"LocalAdaptive"]]},{t,1,21,3}]
ab1=ListPlot[%]
(*ab1=ListPlot[%,{t,0,50},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]*)
Show[p1,ab1]*)
Beep[]

Different attempts (ignore)

Comparison with the swapan gosh for parameters (sink placed at the minima)
g=1; d=.3; x0=5; k0=3; a=0; b=xc; xc=0;

Show[p4,ab2,p1]

(Legended[, "Analytical"])

COmparison with the SK Ghosh answer

Unprotect["Global`*"];
ClearAll["Global`*"]
Remove["Global`*"]
y1[x_]=C1 Exp[(p+q)x];(*x<-b*)
y2[x_]=C3 Exp[(p+q)x]+C4 Exp[(q-p)x];(*0>x>-b*)
y3[x_]=A Exp[(p-q)x]+B*Exp[-(p+q)x];(*x0>x>0*)
(*{{C1,A,B}}={C1,A,B}/.Solve[{y1[0]==y3[0],
de(y3'[0]-y1'[0])+w de y3[0]-w de*y1[0](-1)==-1,de(-y3'[x0]-w y3[x0])==k0*y3[x0]},{C1,A,B}];*)
{{C1,C3,C4,A,B}}={C1,C3,C4,A,B}/.Solve[{y1[-b]==y2[-b],de(y2'[-b]-y1'[-b])==-1,y2[0]==y3[0],de(y3'[0]-y2'[0])+w de y3[0]-w de*y2[0](-1)==0,de(-y3'[x0]-w y3[x0])==k0*y3[x0]},{C1,C3,C4,A,B}]
q=w/2;
p=Sqrt[q^2+(s/de)];
w=2g;(*w=2g*)
de=.3/2;(*de=d/2*)
x0=-1;(*position of the sink*)
k0=3;(*Sink strength*)
Clear[b];
b=4;(*position of the initial condition*)
F[x_,s_]=Piecewise[{{y1[x],x<-b},{y2[x],0> x≥-b},{y3[x],x0≥ x≥0}}];
f[x_,s_]=Integrate[F[x,s],x];
Q[s_]=FullSimplify[f[x0,s]];
(*Print["Check between the answers of survival prob"]
Plot[{ρ[s],Q[s],1/s},{s,1,20},PlotLegends->Automatic]*)
Attributes[NInverseLaplaceTransform]={NumericFunction};

NInverseLaplaceTransform:: wrongexpr="The argument `1` is not Numeric at fixed `2` .";NInverseLaplaceTransform:: wrongtime="The argument `1` is not Numeric.";
NInverseLaplaceTransform:: failed=
"The inverse could not be calculated. Suspect one of the following reasons: \n- the inverse does not exist; \n- the inverse or its derivatives have a discontinuity and the specified time falls near to it; \n- the inverse is zero; \n- the specified time is too large or small when the inverse is periodic-like.";

Options[NInverseLaplaceTransform]={"Startm"->5,"Method"->"Automatic"};

NInverseLaplaceTransform[expr_,s_Symbol,t_?NumericQ, opts:OptionsPattern[]]:=
Module[{F,z, x,M, f,g,fg,gg,ff,fgm, f12,g12, m,mmin,mmax,mb,meth,res }, 

mmin=Min[15,Max[5,OptionValue["Startm"]]];
mmax=Min[20,mmin+5];

meth=OptionValue["Method"];
If[Not[Or[meth==="FT",meth==="GWR"]],meth="Automatic"];

F[z_]=Rationalize[expr,0]/.s->z;
x=Rationalize[t,0];
If[Not[Block[{$MaxExtraPrecision=5000.},NumericQ[F[100]]]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongexpr,expr,s]);res=$Failed;Return[res]];If[Not[And[NumericQ[t],t>0]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongtime,t]);res=$Failed;Return[res]];
f12=g12=Array[0&,mmax];
Do[
M=2^m;mb=m;ff=fgm=fg=gg=False;

If[Or[meth==="FT",meth==="Automatic"],
f=FT[F, x, M]; 
f12[[m]]=If[$MinMachineNumber<Abs[f]<$MaxMachineNumber,f//N,0.];
ff=And[f12[[m]]≠0,f12[[m]]==f12[[m-1]]];
fgm=And[f12[[m]]≠0,f12[[m]]==g12[[m-1]]];
     If[Or[ff,fgm],Break[]];      
     If[And[m==mmax-1,meth==="FT",f12[[m]]==0],res=$Failed;Break[]];
     If[And[m==mmax,f12[[m]]==0,g12[[m-1]]==0],res=$Failed;Break[]];
];

If[Or[meth==="GWR",meth==="Automatic"],
g=GWR[F, x, M];
g12[[m]]=If[$MinMachineNumber<Abs[g]<$MaxMachineNumber,g//N,0.];
fg=And[g12[[m]]≠0,f12[[m]]==g12[[m]]];
gg=And[g12[[m]]≠0,g12[[m]]==g12[[m-1]]];
      If[Or[fg,gg],Break[]];
      If[And[m==mmax-1,meth==="GWR",g12[[m]]==0],res=$Failed;Break[]];   
      If[And[m==mmax,g12[[m]]==0],res=$Failed;Break[]];
];
 If[m==mmax,res=$Failed]
,{m,mmin,mmax}
];
If[Or[fg,gg],res=g12[[mb]]; Return[res]];
If[Or[ff,fgm],res=f12[[mb]]; Return[res]];
res==$Failed;
(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform::failed]);
res
];

GWR[F_, t_?NumericQ, M_:32, precin_:0] := Module[
{ fv,w, wm,wp, best, tau = Log[2]/t,  prec,k,n,i}, 
If[precin ≤ 0, prec = 2.1*M, prec = precin]; 
If[prec ≤ $MachinePrecision, prec = $MachinePrecision]; 
If[Precision[tau] < prec, tau = SetPrecision[tau, prec]]; 
Block[{$MaxExtraPrecision=5000.},
fv=Table[ N[F[i tau], prec], {i, 2M}]; 
w=Table[ 
tau*((2*n)!/(n!*(n - 1)!))*Total[Table[Binomial[n, i]*(-1)^i*fv[[n + i]], {i, 0, n}]]
, {n,  M}]; 
];
best=Last[w];
If[Not[NumberQ[best]],best=$Failed;Return[best]];
wm=Table[0,{M}];
Do[
wp=Quiet[k/Differences[w]];
wp=Quiet[Drop[wm,1]+wp];
wm=Drop[w,-1];
w=wp;
If[EvenQ[k],If[NumberQ[Last[w]],best =Last[w],Break[],Break[]]]
,{k,  M -1}]; 
If[Not[NumberQ[best]],best=$Failed];
best
]

FT[F_,t_?NumericQ,M_:32]:=Module[{np,p,r,s,theta,sigma,tab,f},
np=Max[M,$MachinePrecision];
p=SetPrecision[t,np];
Block[{$MaxExtraPrecision=5000.},
r=2M/(5p);
s=r theta (Cot[theta]+I);
sigma=theta+(theta Cot[theta]-1)Cot[theta];
Quiet[
tab=Flatten[{Exp[p r] F[r]/2,Table[Exp[p s](1+I sigma)F[s],{theta,Pi/M,(M-1)Pi/M,Pi/M}]/.ComplexInfinity->0}];
f=Total[Re[(r/M) tab]]
]
];
f
]
h=-.4;
y=-1.1;
ys=0;
g=0.5
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}];
ab2=ListLinePlot[%,PlotStyle->Dashed,PlotLegends->Placed[{Style["g=0.5",13]},{Scaled[{0,0.5}],{y+1,h-1}}]];
Clear[g]
g=0.8
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}];
ab3=ListLinePlot[%,PlotStyle->{Thick,Black},PlotLegends->Placed[{Style["g=0.8",13]},{Scaled[{0,0.5}],{y+1,h}}]];
Clear[g]
g=1
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}];
ab4=ListLinePlot[%,PlotStyle->Dotted,PlotLegends->Placed[{Style["g=1",13]},{Scaled[{0,0.5}],{y+1,h+1}}],PlotTheme->"Scientific",PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}]
Show[ab4,ab2,ab3]
Beep[]

Let us generate graphs on the effect of g;

ClearAll["Global`*"]
σ=0;
ϵ1=(k1/(2(d/2)));
ϵ2=(k2/(2(d/2)));
k1=-d g;
k2=k0;

P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*UnitStep[b-x]+x0)](-RealSign[((x*RealSign[x-b]+2b*UnitStep[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(x*RealSign[x-b]+2b*UnitStep[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+HeavisideTheta[a-x](((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((x+2(b-a)+x0))/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((-x+2a+x0))/(2Sqrt[(σ^2+ d t)]))]))/.{t->t/2};(*{RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]->1,t->t/2}*)P[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x])-g^2d t/2]*P1[x,t];Qt[t_]=((Exp[g(RealAbs[x0]-x0)])/2)(Erfc[((Sqrt[(-x0+d g t)^2])/(Sqrt[2 d t]))]+Exp[2 g x0]Erfc[((d g t+x0)/(Sqrt[2 d t]))])+((Exp[ϵ2^2d t/2+ϵ2 x0+g RealAbs[x0]-g^2 d t/2])/2)(ϵ2 Exp[ϵ2*2b](-RealSign[((2b+x0)/(d t))+ϵ2]+Erf[((2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])((1/(g-ϵ2))+(1/(g+ϵ2)))+ϵ2 Exp[(ϵ2-g)b](-RealSign[((b+x0)/(d t))+ϵ2]+Erf[((b+x0+ϵ2 d t)/(Sqrt[2 d t]))])(-(1/(ϵ2-g))-(1/(ϵ2+g)))-ϵ2 Exp[ϵ2(-x0-ϵ2 d t)]((2/(g-ϵ2))Exp[g(2b+x0+ϵ2 d t)]UnitStep[-2b-x0-ϵ2 d t]+(2/(-(ϵ2+g)))Exp[-g(2b+x0+ϵ2 d t)]UnitStep[2b+x0+ϵ2 d t]UnitStep[b-2b-x0-ϵ2 d t]-((2 Exp[g(x0+ϵ2 d t)])/(ϵ2-g))UnitStep[-x0-ϵ2 d t-b])+ϵ2(((1/(g-ϵ2)))Exp[2b g+((g-ϵ2)/2)(2x0+(g+ϵ2)d t)]*Erfc[((2b+x0+d g t)/(Sqrt[2 d t]))]-((Exp[((g-ϵ2)/2)(2x0+(g+ϵ2)d t)])/(ϵ2-g))*Erfc[((b+d g t+x0)/(Sqrt[2 d t]))]+((Exp[-2b g+(1/2)(-2x0+d t(g-ϵ2))(g+ϵ2)])/(g+ϵ2))(Erfc[((Sqrt[(2b-d g t+x0)^2])/(Sqrt[2 d t]))]-Erfc[((Sqrt[(b-d g t+x0)^2])/(Sqrt[2 d t]))]))+ϵ1*Exp[2b ϵ2](((Exp[-ϵ2 a+g a])/(g-ϵ2))(-Sign[((-a+2b+x0)/(d t))+ϵ2]+Erf[((-a+2b+x0+ϵ2 d t)/(Sqrt[2 d t]))])-(2/(g-ϵ2))Exp[(g-ϵ2)(2b+x0+ϵ2 d t)]UnitStep[a-2b-x0-ϵ2 d t]+(2/(2(g-ϵ2)))Exp[((g-ϵ2)/2)(4b+2x0+d t(g+ϵ2))]*Erfc[((2b-a+d g t+x0)/(Sqrt[2d t]))])-ϵ1*Exp[ϵ2 2(b-a)](((Exp[(ϵ2+g) a])/(g+ϵ2))(-Sign[((a+2(b-a)+x0)/(d t))+ϵ2]+Erf[((a+2(b-a)+x0+ϵ2 d t)/(Sqrt[2 d t]))])-((Exp[-(g+ϵ2)(2(b-a)+x0+ϵ2 d t)])/(g+ϵ2))(-2)UnitStep[-a+2b+x0+ϵ2 d t]-((Exp[((g+ϵ2)/2)(4a-4b+d g t-x0-ϵ2 d t)])/(g+ϵ2))Erfc[(((a-2b+d g t-x0))/(Sqrt[2d t]))]));g=.1;
d=.3;
x0=0;
k0=3;
a=0;
b=xc;
xc=-1;
h=-.4;
y=-1.1;
ys=0;
Ex[t_]:=0ϵ1/2*NIntegrate[Erfc[((Sqrt[(x+x0)^2])/(Sqrt[2 d t]))]+Erfc[((Sqrt[(-x+2b+x0)^2])/(Sqrt[2d t]))]-Erfc[((Sqrt[(x+2(b-a)+x0)^2])/(Sqrt[2d t]))]-Erfc[((Sqrt[(-x+2a+x0)^2])/(Sqrt[2d t]))],{x,-∞,a},Method->"AdaptiveMonteCarlo",MaxRecursion->20]
(*c=f/.Flatten[Block[{t=0.1},Solve[f*(Qt[t]+Ex[t])==1,f]]];
Table[c(Qt[t]+Ex[t]),{t,0.1,50}]
p1=ListLinePlot[%,PlotRange->Full,PlotStyle->Gray,PlotLegends->Placed[{Style["g=0.1",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotTheme->{"Scientific"}];*)
ClearAll[c,g]
g=.5;
c=f/.Flatten[Block[{t=0.1},Solve[f*(Qt[t]+Ex[t])==1,f]]];
Table[c(Qt[t]+Ex[t]),{t,0.1,50}]
p2=ListLinePlot[%,PlotRange->Full,PlotStyle->Dashed,PlotLegends->Placed[{Style["g=0.5",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotTheme->"Scientific"];
Table[c NIntegrate[P1[x,t],{x,-30,30}],{t,0.1,50}];
ab1=ListLinePlot[%]
ClearAll[c,g]
g=.8;
c=f/.Flatten[Block[{t=0.1},Solve[f*(Qt[t]+Ex[t])==1,f]]];
Table[c(Qt[t]+Ex[t]),{t,0.1,50}]
p3=ListLinePlot[%,PlotRange->Full,PlotStyle->{Thick},PlotLegends->Placed[{Style["g=0.8",13]},{Scaled[{0,0.5}],{y-4.5ys,h}}],PlotTheme->{"Scientific","Monochrome"}];
Table[c NIntegrate[P1[x,t],{x,-30,30}],{t,0.1,50}];
ab2=ListLinePlot[%]
ClearAll[c,g]
g=1;
c=f/.Flatten[Block[{t=0.1},Solve[f*(Qt[t]+Ex[t])==1,f]]];
Table[c(Qt[t]+Ex[t]),{t,0.1,50}]
p4=ListLinePlot[%,PlotRange->Full,PlotStyle->{Dotted,Black},PlotLegends->Placed[{Style["g=1",13]},{Scaled[{0,0.5}],{y-4.5ys,h+1}}],PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]}];
Table[c NIntegrate[P1[x,t],{x,-30,30}],{t,0.1,50}];
ab3=ListLinePlot[%]
Show[p4,p2,p3]
Show[ab1,ab2,ab3]

ClearAll["Global`*"];
g=.7;
d=.3;
x0=3;
k0=3;
a=0;
b=xc;
xc=1;
h=-.4;
y=-1.1;
ys=0;
Ex[t_]:=NIntegrate[Erfc[((x+x0)/(Sqrt[2 d t]))]+Erfc[((-x+2b+x0)/(Sqrt[2d t]))]-Erfc[((x+2(b-a)+x0)/(Sqrt[2d t]))]-Erfc[((-x+2a+x0)/(Sqrt[2d t]))],{x,-∞,a}]
Plot[Ex[t],{t,0,50}]

Show[ab4,ab2,ab3]

When sink is at 0, it matches with sink placed at 0 solution;

Following similar calculation strategies as done in the absolute potential version notebook.
Solution using our transformations for the problem without the sink by random variable transformation

ClearAll["Global`*"]
a=1;
b=0;
Plot[2(UnitStep[b-x]RealAbs[x+a]-UnitStep[x-b](x-a)),{x,-3,3},PlotRange->Full]

ClearAll["Global`*"];
σ=0;
(*FX[x1_,t_]=Assuming[x≠0,Integrate[(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2,{x,-∞,x1}]];*)
P[x_,t_]=(Exp[-g*(UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)-(UnitStep[x0]RealAbs[-x0+a]+UnitStep[-x0](a+x0)))-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;(*D[FX[x1,t],x1]*)
Print["Analytical check for the total solution"]
Assuming[x≠-a&&x≠0&&x<-a,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*(UnitStep[-x]RealSign[x+a]-UnitStep[x]),x]-D[P[x,t],t]]]
Assuming[x≠-a&&-a<x<a&&x≠0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*(UnitStep[-x]RealSign[x+a]-UnitStep[x]),x]-D[P[x,t],t]]]
Assuming[x≠-a&&x>a&&x≠0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t](UnitStep[-x]RealSign[x+a]-UnitStep[x]),x]-D[P[x,t],t]]]
Pmaplus[x_,t_]=P[x,t]/.UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)->x+a;
Pmaminus[x_,t_]=P[x,t]/.UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)-> -x-a;
P0plus[x_,t_]=P[x,t]/.UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)-> a-x;
P0minus[x_,t_]=P[x,t]/.UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)->x+a;
Print["Intial condition check"]
P[-x0,t]
Print["Continuity check"];
Simplify[Pmaplus[-a,t]-Pmaminus[-a,t]]
Print["Boundary check of Mr. R Saravanan which leads to differential equation at x=-a"]
Simplify[d/2 D[Pmaplus[x,t],{x,1}]+d g*Pmaplus[x,t]-d/2 D[Pmaminus[x,t],{x,1}]-d g*Pmaminus[x,t]*(-1)-d g Pmaplus[-a,t]/.x->-a]
Print["Continuity check"];
Simplify[P0plus[0,t]-P0minus[0,t]]
Print["Boundary check of Mr. R Saravanan which leads to differential equation at x=0"]
Simplify[d/2 D[P0plus[x,t],{x,1}]+d g*P0plus[x,t](-1)-d/2 D[P0minus[x,t],{x,1}]-d g*P0minus[x,t]+d  g P0plus[0,t]/.x->0]

ClearAll["Global`*"];
σ=0;P[x_,t_]=(Exp[-g*(UnitStep[-x]RealAbs[x+a]+UnitStep[x](a-x)-(UnitStep[x0]RealAbs[-x0+a]+UnitStep[-x0](a+x0)))-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
x0=10;
g=.3;
d=0.5;
a=-2;
Plot[{P[x,1],P[x,10],P[x,30],P[x,100]},{x,-20,20},PlotRange->Full](*D[FX[x1,t],x1]*)

But as usual it is not going to be normalized/follow all the boundary conditions; Hence we add the other boundary effects through following equation

-d g (DiracDeltaSeq[x]) P(x,t)+d g (DiracDeltaSeq[x+a])P(x,t)+(1/2) d Derivative[2,0][P](x,t)==Derivative[0,1][P](x,t)

THe model is similar to that of the Oster-Nishijima’s model. And the solution is already derived by us. Let us do all the tests before analyzing numerically.

Note that no term should be discarded if you want all the checks to be satisfied

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
Print["Our solution is"]
σ=0;
ϵ1=(k1/(2(d/2)));
ϵ2=(k2/(2(d/2)));
k1=-d g;
k2=d g;
Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
Heaviside[x_]:=Piecewise[{{0,x≤ 0},{1,x>0}}];
P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](-RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+Heaviside[a-x](((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((x+2(b-a)+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((-x+2a+x0)/(2Sqrt[(σ^2+ d t)]))]))/.{RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]->1,t->t/2};Pe[x_,t_]=Exp[-g*(Unit[-x]RealAbs[x+xm]+Unit[x](xm-x)-(UnitStep[x0]RealAbs[-x0+xm]+UnitStep[-x0](xm+x0)))-g^2d* t/2]*P1[x,t];
Pebplus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]->x,Heaviside[a-x]->0,Unit[-x]RealAbs[x+xm]+Unit[x](xm-x)->xm-x};
Pebminus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->0,Unit[-x]RealAbs[x+xm]+Unit[x](xm-x)->x+xm};
Peaplus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->0,Unit[-x]RealAbs[x+xm]+Unit[x](xm-x)->x+xm};
Peaminus[x_,t_]=c Pe[x,t]/.{x*RealSign[x-b]+2b*Unit[b-x]-> 2b-x,Heaviside[a-x]->1,Unit[-x]RealAbs[x+xm]+Unit[x](xm-x)-> -x-xm};
a=-xm;
b=0;
Print["Solubility check of Pe"]
Assuming[x≠ a&&x≠b&&x0>0&&xm>0,Simplify[d/2 D[Pe[x,t],{x,2}]+d g* D[Pe[x,t]*(Unit[-x]RealSign[x+xm]-Unit[x]),x]-D[Pe[x,t],t]]]
Print["Continuity check of P3 a, b"]
Simplify[Pebplus[b,t]-Pebminus[b,t]]
Simplify[Peaplus[a,t]-Peaminus[a,t]]
Print["Boundary check for Pe"]
Assuming[b≥a,Simplify[d/2 D[Peaplus[x,t],{x,1}]+d g*Peaplus[x,t]-d/2 D[Peaminus[x,t],{x,1}]-d g*Peaminus[x,t]*(-1)/.x->a]]
Assuming[b≥a,Simplify[d/2 D[Pebplus[x,t],{x,1}]-d g*Pebplus[x,t]-d/2 D[Pebminus[x,t],{x,1}]-d g*Pebminus[x,t]/.x->b]]
Beep[]

Constructing the solution for same problem using Privman’s method

ClearAll["Global`*"];
y1[x_]=A*Exp[(q+r)x/a]
y1b[x_]=A1*Exp[(q+r)x/a]+b1*Exp[-q x/a]
y2[x_]=c1*Exp[(q)x/a]+c2*Exp[-(q+r)x/a]
y3[x_]=B*Exp[-(q)x/a]
{{A,A1,b1,c1,c2,B}}={A,A1,b1,c1,c2,B}/.Solve[{y1[-b]==y1b[-b],y1b[-a]==y2[-a],y2[0]==y3[0],de(y2'[-a]-y1b'[-a])+2r de/a y2[-a]==0,de(y1b'[-b]-y1'[-b])==1,de(y3'[0]-y2'[0])-2r de/a y2[0]==0},{A,A1,b1,c1,c2,B}]
z=((4a^2)/(r^2 de))s;
q=(r/2)(Sqrt[1+z]-1);
F[x_,s_]=Piecewise[{{y1[x],x<-b},{y1b[x],-a>x≥ -b},{y2[x],0>x≥ -a},{y3[x],x≥ 0}}];
de=.4;
b=4;
a=3;
g=0.3;
r=2a g;
Print["Normalisation check"]
f[x_,s_]=Integrate[F[x,s],x];
Simplify[f[x,s]/.x->-∞];
Simplify[f[x,s]/.x->∞];
Normm=Simplify[%%-%]
Q[s_]=-Simplify[f[0,s]/s,Assumptions->b>a&&b∈Positive];

We have got to check all the tests, Now numerical test by taking inverse Laplace transform of V Privman solutions. THe comparison with the SK gosh solution of Privman is tested below;

Unprotect["Global`*"];
ClearAll["Global`*"]
Remove["Global`*"]
Print["Our solution is"]
σ=0;
Unit[x_]:=Piecewise[{{0,x<0},{1,x>0},{.5,x==0}}];
Heaviside[x_]:=Piecewise[{{0,x≤ 0},{1,x>0}}];
P1[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](-RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[(x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+Heaviside[a-x](((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((x+2(b-a)+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((-x+2a+x0)/(2Sqrt[(σ^2+ d t)]))]))/.{RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x*RealSign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+ϵ2]->1,RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+ϵ2]->1,t->t/2};
Pnc[x_,t_]:=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]-((ϵ2)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*RealSign[x-b]+2b*Unit[b-x]+x0)](Erfc[(x*RealSign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)])])+Heaviside[a-x](((-ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](Erfc[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](Erfc[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((x+2(b-a)+x0))/(2Sqrt[(σ^2+ d t)]))]-Erfc[(((-x+2a+x0))/(2Sqrt[(σ^2+ d t)]))]))/.{t->t/2};P[x_,t_]=Exp[-g*(UnitStep[-x]RealAbs[x+xm]+UnitStep[x](xm-x)-(UnitStep[x0]RealAbs[-x0+xm]+UnitStep[-x0](xm+x0)))-g^2d* t/2]*Pnc[x,t];
ϵ1=(k1/(2(d/2)));
ϵ2=(k2/(2(d/2)));
k1=-d g;
k2=d g;
a=-xm;
xm=2;(*potential minimum in positive number*)
b=0;(*PLease dont change, this is assumed for the privman solution*)
g=.3;
d=.3;
x0=5;(*position of the initial condition*)
h=-.4;
y=-1.1;
ys=0;
f=c/.Solve[c*NIntegrate[P[x,.1],{x,-30,30},AccuracyGoal->Infinity,MaxRecursion->30,Method->"AdaptiveQuasiMonteCarlo"]==1,c]
p2=Plot[{f P[x,.1],f P[x,1],f P[x,10],f P[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic]
Flatten[Table[f NIntegrate[ P[x,t],{x,-30,0},AccuracyGoal->Infinity,MaxRecursion->30,Method->"AdaptiveQuasiMonteCarlo"],{t,.01,50}]]
p1=ListPlot[%,PlotJoined->True,PlotRange->Full,PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}]]
Attributes[NInverseLaplaceTransform]={NumericFunction};

NInverseLaplaceTransform:: wrongexpr="The argument `1` is not Numeric at fixed `2` .";NInverseLaplaceTransform:: wrongtime="The argument `1` is not Numeric.";
NInverseLaplaceTransform:: failed=
"The inverse could not be calculated. Suspect one of the following reasons: \n- the inverse does not exist; \n- the inverse or its derivatives have a discontinuity and the specified time falls near to it; \n- the inverse is zero; \n- the specified time is too large or small when the inverse is periodic-like.";

Options[NInverseLaplaceTransform]={"Startm"->5,"Method"->"Automatic"};

NInverseLaplaceTransform[expr_,s_Symbol,t_?NumericQ, opts:OptionsPattern[]]:=
Module[{F,z, x,M, f,g,fg,gg,ff,fgm, f12,g12, m,mmin,mmax,mb,meth,res }, 

mmin=Min[15,Max[5,OptionValue["Startm"]]];
mmax=Min[20,mmin+5];

meth=OptionValue["Method"];
If[Not[Or[meth==="FT",meth==="GWR"]],meth="Automatic"];

F[z_]=Rationalize[expr,0]/.s->z;
x=Rationalize[t,0];
If[Not[Block[{$MaxExtraPrecision=5000.},NumericQ[F[100]]]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongexpr,expr,s]);res=$Failed;Return[res]];If[Not[And[NumericQ[t],t>0]],(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform:: wrongtime,t]);res=$Failed;Return[res]];
f12=g12=Array[0&,mmax];
Do[
M=2^m;mb=m;ff=fgm=fg=gg=False;

If[Or[meth==="FT",meth==="Automatic"],
f=FT[F, x, M]; 
f12[[m]]=If[$MinMachineNumber<Abs[f]<$MaxMachineNumber,f//N,0.];
ff=And[f12[[m]]≠0,f12[[m]]==f12[[m-1]]];
fgm=And[f12[[m]]≠0,f12[[m]]==g12[[m-1]]];
     If[Or[ff,fgm],Break[]];      
     If[And[m==mmax-1,meth==="FT",f12[[m]]==0],res=$Failed;Break[]];
     If[And[m==mmax,f12[[m]]==0,g12[[m-1]]==0],res=$Failed;Break[]];
];

If[Or[meth==="GWR",meth==="Automatic"],
g=GWR[F, x, M];
g12[[m]]=If[$MinMachineNumber<Abs[g]<$MaxMachineNumber,g//N,0.];
fg=And[g12[[m]]≠0,f12[[m]]==g12[[m]]];
gg=And[g12[[m]]≠0,g12[[m]]==g12[[m-1]]];
      If[Or[fg,gg],Break[]];
      If[And[m==mmax-1,meth==="GWR",g12[[m]]==0],res=$Failed;Break[]];   
      If[And[m==mmax,g12[[m]]==0],res=$Failed;Break[]];
];
 If[m==mmax,res=$Failed]
,{m,mmin,mmax}
];
If[Or[fg,gg],res=g12[[mb]]; Return[res]];
If[Or[ff,fgm],res=f12[[mb]]; Return[res]];
res==$Failed;
(ResourceFunction["ResourceFunctionMessage"][NInverseLaplaceTransform::failed]);
res
];

GWR[F_, t_?NumericQ, M_:32, precin_:0] := Module[
{ fv,w, wm,wp, best, tau = Log[2]/t,  prec,k,n,i}, 
If[precin ≤ 0, prec = 2.1*M, prec = precin]; 
If[prec ≤ $MachinePrecision, prec = $MachinePrecision]; 
If[Precision[tau] < prec, tau = SetPrecision[tau, prec]]; 
Block[{$MaxExtraPrecision=5000.},
fv=Table[ N[F[i tau], prec], {i, 2M}]; 
w=Table[ 
tau*((2*n)!/(n!*(n - 1)!))*Total[Table[Binomial[n, i]*(-1)^i*fv[[n + i]], {i, 0, n}]]
, {n,  M}]; 
];
best=Last[w];
If[Not[NumberQ[best]],best=$Failed;Return[best]];
wm=Table[0,{M}];
Do[
wp=Quiet[k/Differences[w]];
wp=Quiet[Drop[wm,1]+wp];
wm=Drop[w,-1];
w=wp;
If[EvenQ[k],If[NumberQ[Last[w]],best =Last[w],Break[],Break[]]]
,{k,  M -1}]; 
If[Not[NumberQ[best]],best=$Failed];
best
]

FT[F_,t_?NumericQ,M_:32]:=Module[{np,p,r,s,theta,sigma,tab,f},
np=Max[M,$MachinePrecision];
p=SetPrecision[t,np];
Block[{$MaxExtraPrecision=5000.},
r=2M/(5p);
s=r theta (Cot[theta]+I);
sigma=theta+(theta Cot[theta]-1)Cot[theta];
Quiet[
tab=Flatten[{Exp[p r] F[r]/2,Table[Exp[p s](1+I sigma)F[s],{theta,Pi/M,(M-1)Pi/M,Pi/M}]/.ComplexInfinity->0}];
f=Total[Re[(r/M) tab]]
]
];
f
]
ClearAll[z,q,Ω,de,a,r,b,F,f,s]
y1[x_]=A*Exp[(q+r)x/a](*x<-b*)
y1b[x_]=A1*Exp[(q+r)x/a]+b1*Exp[-q x/a](*x<-a*)
y2[x_]=c1*Exp[(q)x/a]+c2*Exp[-(q+r)x/a](*-a<x<0*)
y3[x_]=B*Exp[-(q)x/a](*x>0*)
{{A,A1,b1,c1,c2,B}}={A,A1,b1,c1,c2,B}/.Solve[{y1[-b]==y1b[-b],y1b[-a]==y2[-a],y2[0]==y3[0],de(y2'[-a]-y1b'[-a])+2r de/a y2[-a]==0,de(y1b'[-b]-y1'[-b])==-1,de(y3'[0]-y2'[0])-2r de/a y2[0]==0},{A,A1,b1,c1,c2,B}]
z=((4a^2)/(r^2 de))s;
q=(r/2)(Sqrt[1+z]-1);
F[x_,s_]=Piecewise[{{y1[x],x<-b},{y1b[x],-a>x≥ -b},{y2[x],0>x≥ -a},{y3[x],x≥ 0}}];
de=d/2;
a=2;(*position of the minima given in positive number*)
b=x0;(*position of initial condition*)
r=2a g;
Print["Normalisation check"]
f[x_,s_]=Integrate[F[x,s],x];
Simplify[f[x,s]/.x->-∞];
Simplify[f[x,s]/.x->∞];
Normm=Simplify[%%-%]
Q[s_]=Simplify[f[0,s]-f[-30,s],Assumptions->b≥ a&&b∈Positive]
Ω=4q(r+q)Exp[r+2q]+r^2;
ρ[s_]=((a^2)/((r+q)de Ω))(4(r+q)Exp[r+2q]+r^2*((1-Exp[q])/q)-(3r+2q)Exp[q]);
Table[NInverseLaplaceTransform[Q[s],s,t],{t,.1,50}]
ab2=ListPlot[%,PlotLegends->"Privmanb",PlotStyle->Green]
Table[NInverseLaplaceTransform[ρ[s],s,t],{t,.1,50}]
ab1=ListPlot[%,PlotLegends->"Privmanhim",PlotStyle->Red]
Show[p1,ab1,ab2]
Beep[]

Show[ab2,p1]

(Legended[, ({{"Privmanhim"}, {"Privmanb"}})])

Thus we derived a new method of converting piecewise linear potential with two-boundary problems as an Oster-Nishijima model. Using the method, we derived an exaact solution of P(x,t). The survival probability is exact when the sink is placed at the minima. The disadvantage ofthe solution is the artifact of Oster-Nishijima solution.

Miscellaneous trials.

ClearAll["Global`*"];
g=1;
d=.3;
x0=5;
k0=3;
a=0;
b=xc;
xc=-2;
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
dThick=20;
dTime=50;
eqSimplee=D[Pe[x,t],t]==d/2*D[Pe[x,t],{x,2}]+d g RealSign[x]*D[Pe[x,t],x]+2d g S[x]*Pe[x,t]-k0*S[x-xc]*Pe[x,t];
eqInie=Pe[x,0]==√(1000/π)*Exp[-1000*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
h=-.4;
y=-1.5;
ys=0;
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqInie,eqBlefte,eqBleftg},{Pe[x,t]},{x,-dThick,dThick},{t,0,dTime},MaxStepSize->0.001]
Flatten[Table[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick},AccuracyGoal->10,PrecisionGoal->6,Method->"LocalAdaptive"],{t,1,50}]]
ab1=ListPlot[%]
(*ab1=ListPlot[%,{t,0,50},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]*)
Show[ab1]

Table[{t,t^2},{t,0,100}]

ListPlot[%]

D[RealSign[x],x]

Integrate[RealSign[x],{x,-2,2}]
Integrate[2UnitStep[x]-1,{x,-2,2}]
Integrate[D[RealSign[x],x],{x,-2,2}]
Integrate[D[2UnitStep[x]-1,x],{x,-2,2}]
Integrate[2DiracDelta[x],{x,-2,2}]

D[2UnitStep[x]-1,x]

DiracDelta[x]

Integrate[UnitStep'[x],{x,-2,2}]

ClearAll["Global`*"];
Γ= g d;
l=d/(2Γ);
σ=0;
u[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))];
Print["Analytical check of kenkre"]
(*Assuming[x≠a,Simplify[d/2 D[u[x,t],{x,2}]+d g* D[u[x,t]*RealSign[x-a],x]-D[u[x,t],t]]]*)
uplus[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))]/.
{x*RealSign[x]->x};
uminus[x_,t_]=(Exp[-1/(2l)(x*RealSign[x]-x0 RealSign[x0])]/Sqrt[4π(σ^2+d/2 t)])Exp[((-(x+x0)^2-Γ^2 t^2)/(4(σ^2+d/2 t)))]+((Exp[-((x*RealSign[x])/l)])/(4l))*Erfc[((x*RealSign[x]+x0 RealSign[x0]-Γ t)/(Sqrt[4 (σ^2+d/2t)]))]/.{x*RealSign[x]->-x};
Print["Continuity check"]
Simplify[Limit[uplus[x,t]-uminus[x,t],x->0]]
Print["Initial condition check"]
x=-x0;
u[x,t]
Clear[x];
Print["Boundary check of kenkre"]
Simplify[Limit[d/2 D[uplus[x,t],{x,1}]+d g*uplus[x,t]-d/2 D[uminus[x,t],{x,1}]-d g*uminus[x,t]*(-1),x->0]]
Simplify[d/2 D[uplus[x,t],{x,1}]+d g*uplus[x,t]-d/2 D[uminus[x,t],{x,1}]-d g*uminus[x,t]*(-1)/.x->0]
P[x_,t_]=(Exp[-g*((x*RealSign[x])-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Pplus[x_,t_]=(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.{t->t/2,(x*RealSign[x])->x};
Pminus[x_,t_]=(Exp[-g*(x*RealSign[x]-x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.{t->t/2,(x*RealSign[x])->-x};
P3[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(x*RealSign[x]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((x*RealSign[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((x*RealSign[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]))/.RealSign[((x*RealSign[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1;
Simplify[P[x,t]]
Simplify[P3[x,t]]
Print["Analytical check of Mr. R Saravanan"]
Assuming[x≠ 0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
σ=0;
Print["Continuity check"]
Simplify[Pplus[0,t]-Pminus[0,t]]
Print["Boundary check of Mr. R Saravanan which leads to differential equation at x=0"]
Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)-d g Pplus[0,t]/.x->0]
(*Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)/.x->a]
Simplify[-d g Pplus[a,t]/.x->a]
Simplify[-d g P[a,t]/.x->a]*)
Print["The analytical difference is"]
Simplify[u[x,t]-P[x,t],Assumptions->Element["Global`*",Reals]&&Element[g,Reals]&&g>0&&Element[d,Reals]&&d>0&&Element[t,Reals]&&t>0]
g=1;
d=.1;
x0=5;
a=-2;
Plot[{u[x,50],P[x,50]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic]
p1=Plot[{u[x,.1],u[x,1],u[x,10],u[x,500]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic]
p2=Plot[{P[x,.1],P[x,1],P[x,10],P[x,500]},{x,-20,20},PlotRange->Full,PlotLegends->Automatic]
Show[p1,p2]

TraditionalForm[ToExpression["\\lim_{\\epsilon\\to0}",TeXForm]]
TraditionalForm[ToExpression["\\frac{D}{2}\\frac{\\partial P(x,t)}{\\partial x}+gDP(x,t)\\text{sgn}(x)",TeXForm]]
TraditionalForm[ToExpression["D g P[0,t]",TeXForm]]

D g sgn(x) P(x,t)+(1/2) D Derivative[1,0][P](x,t)//StandardForm

ClearAll["Global`*"];
Pp[x_,t_]=(Exp[-g(x)+g(x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Print["Analytical check for x>0"]
Simplify[d/2 D[Pp[x,t],{x,2}]+d g D[Pp[x,t],x]-D[Pp[x,t],t]]
Pm[x_,t_]=(Exp[-g(-x+2a)+g(x0*RealSign[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Print["Analytical check for x<0"]
Simplify[d/2 D[Pm[x,t],{x,2}]-g d*D[Pm[x,t],x]-(∂_t)Pm[x,t]]
Pn[x_,t_]=(Exp[0g(a+ (x0*RealSign[x0]))-0g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Print["Analytical check for x==0"]
Simplify[d/2 D[Pn[x,t],{x,2}]+0g d*D[Pp[x,t],x]-(∂_t)Pn[x,t]]
P[x_,t_]=(Exp[-g*(x*RealSign[x]-(x0*RealSign[x0]))-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
Print["Analytical check of the total solution in piecewise regions x>0 and x<0"](*However the results doesnt depend on the values of x0 and a*)
Assuming[x≠0&&x>0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
Assuming[x≠0&&x< 0,Simplify[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t]]]
Limit[d/2 D[P[x,t],{x,2}]+d g* D[P[x,t]*RealSign[x],x]-D[P[x,t],t],x->0]
Print["Initial check"]
x=-x0;
Simplify[P[x,0]]
Clear[x]
Print["Checking the normalisation"]
Simplify[Integrate[P[x,t],{x,-∞,∞}],Assumptions->Element[σ,Positive]&&Element[g,Positive]&&Element[d,Positive]&&Element[t,Positive]&&Element[a,Reals]&&x0>0&&Re[d t+2 (σ^2)]>0]
Beep[]
Print["If it is not normalised it is due to half part of the solution"]

ToExpression["
P_0(x,t)=\\frac{e^{\\left(-g (x\\sign{x}-x_0\\sign{x_0})-D g^2\\frac{t}{2}\\right)}e^{-\\frac{(x+x_0)^2}{4 D\\frac{t}{2}}}}{2\\sqrt{\\pi D\\frac{t}{2}}}",TeXForm]
TraditionalForm[%]

ToExpression["\\tilde{P}(\\tilde{x},\\tilde{t})=v(x,t)P_0(x(\\tilde{x}),t(\\tilde{t}))",TeXForm]
ToExpression["\\tilde{x}=x \\\\\\tilde{t}=\\frac{t}{2}\\\\
v(x,t)=e^{\\left(-g (x\\text{sgn}(x))-D g^2 \\frac{t}{2}\\right)}",TeXForm]
TraditionalForm[%%]
TraditionalForm[%%]

the Eq.3 becomes,

ClearAll["Global`*"];TraditionalForm[ToExpression["\\frac{D}{2}\\frac{\\partial^2 P_0(x,t)}{\\partial x^2}+gD\\frac{\\partial P_0(x,t)}{\\partial x}=\\frac{\\partial P_0(x,t)}{\\partial t}",TeXForm]]
TraditionalForm[ToExpression["\\frac{D}{2}\\frac{\\partial^2 P_0(x,t)}{\\partial x^2}-gD\\frac{\\partial P_0(x,t)}{\\partial x}=\\frac{\\partial P_0(x,t)}{\\partial t}",TeXForm]]

Before correction it was like this

(Legended[, ({{(LineLegend[1, 2, 3])}, {(LineLegend[1, 2, 3])}})])

1	Samanta, A. and Ghosh, S.K., 1992. Exact results on diffusion from a piecewise linear potential well. The Journal of chemical physics, 97(12), pp.9321-9323.
	2	Privman, V. ed., 2005. Nonequilibrium statistical mechanics in one dimension. Cambridge University Press.
	3	Privman, V. and Frisch, H.L., 1991. Exact solution of the diffusion problem for a piecewise linear barrier. The Journal of chemical physics, 94(12), pp.8216-8219.

-(Subscript[k, 0])S(x)P(x,t)

((∂P (x, t))/(∂t))  =D(((∂^2)P)/(∂(x^2)))+(1/((Subscript[k, B])T)) (∂/(∂x))(P((∂U)/(∂x))),---(1)

Solution using our transformations for the problem without the sink by random variable transformation

Solution using our transformations for the problem without the sink by mere substitution

ClearAll[x,a];
(x*RealSign[x-a]+2a*UnitStep[-x+a])/.x->a
Limit[(x*RealSign[x-a]+2a*UnitStep[-x+a]),x->a]
(*2UnitStep[0.]-1
RealSign[0]
Sign[0]
HeavisideTheta[0]
RealSign[0]*)

ClearAll["Global`*"]
y1[x_]=a1*Exp[I k x]+b Exp[-I k x];
y2[x_]=c Exp[I k x];
p0[x_]=Exp[I k x]+0(((-I g)/(k+I g)))Exp[I k Abs[x]];
a1=1;
{{b,c}}={b,c}/.Simplify[Solve[{y1[a]==y2[a],(d/2)(y2'[a]-y1'[a])==-2d g*y2[a]+d g*p0[a]},{b,c}]]

Clear[P,d,g]
(d/2)D[P[x,t],{x,2}]+d g P[x,t]DiracDelta[x]==D[P[x,t],t]//TraditionalForm

We found the exact solution for sink placed at minimum of the potential;Let us analyse it

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]

Constructing the solution with sink at minimum of the potential

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
Print["Our solution is"]
σ=0;
ϵ0=((-d g)/(2(d/2)));
P1[x_,t_]:=(Exp[-g*(x*RealSign[x-a]+2a*UnitStep[-x+a]-RealAbs[x0])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;P2[x_,t_]=((ϵ0)/2)Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])-g^2d t/2]*Exp[ϵ0(x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])+ϵ0^2d/2t](-RealSign[(1/(2d t/2))((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0])+ϵ0]+Erf[((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0]+2ϵ0 d/2 t)/(2Sqrt[d t/2])])/.RealSign[(1/(2d t/2))((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0])+ϵ0]->1;
(*,Assumptions->Element[d,Positive]&&Element[g,Positive]&&Element[t,Positive]&&Element[x0,Positive]];*)
b=(( k0)/(2(d/2)));
P3[x_,t_]=(b/2)*Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])-g^2d t/2]*Exp[b(RealAbs[x]+RealAbs[x0])+b^2(σ^2+d/2 t)]*(-RealSign[((RealAbs[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((RealAbs[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))])/.RealSign[((RealAbs[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1;
P3aplus[x_,t_]=P3[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->x};
P3aminus[x_,t_]=P3[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->2a-x};
P3plus[x_,t_]=P3[x,t]/.{RealAbs[x]->x};
P3minus[x_,t_]=P3[x,t]/.{RealAbs[x]->-x};
Print["Solubility check of P3"]
Assuming[x≠ a,Simplify[d/2 D[P3[x,t],{x,2}]+d g* D[P3[x,t]*RealSign[x-a],x]-D[P3[x,t],t]]]
Simplify[P3plus[0,t]-P3minus[0,t]]
Print["Continuity check of P3 a, 0"]
Simplify[P3aplus[a,t]-P3aminus[a,t]]
Simplify[P3plus[0,t]-P3minus[0,t]]
Print["Boundary check for P3"]
Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&((k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[P3plus[x,t],{x,1}]+d g*P3plus[x,t]-d/2 D[P3minus[x,t],{x,1}]-d g*P3minus[x,t]*(-1)-( 2d  g +k0)P3minus[0,t]-k0 P1[0,t]/.x->0]]
Print["Our solution is"]
P[x_,t_]=P1[x,t]+P2[x,t]+P3[x,t]
(*Here we start assuming a to be negative*)
Pplus[x_,t_]=P[x,t]/.{RealAbs[x]->x};
Pminus[x_,t_]=P[x,t]/.{RealAbs[x]->-x};
Paplus[x_,t_]=P[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->x};
Paminus[x_,t_]=P[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->2a-x};
Print["Continuity check for the total solution at a, 0"]
Simplify[Paplus[a,t]-Paminus[a,t]]
Simplify[Pplus[0,t]-Pminus[0,t]]
Print["Boundary check at x=a leads to differential equation"]
Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&(( k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[Paplus[x,t],{x,1}]+d g*Paplus[x,t]-d/2 D[Paminus[x,t],{x,1}]-d g*Paminus[x,t]*(-1)-g d P3aminus[a,t]/.x->a]]
Print["Boundary at x=0 leads to differential equation"]
Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&(( k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)-(2d g + k0)Pplus[0,t]+k0*P2[0,t]/.x->0]]
Beep[]

After modelling the boundary conditions as a PDE. and we add corrections to findP4 and P5; lets see what happens; THe modification works

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
Print["Our solution is"]
σ=0;
ϵ0=((-d g)/(2(d/2)));
P1[x_,t_]:=(((Exp[-g*(x*RealSign[x-a]+2a*UnitStep[-x+a]-RealAbs[x0])-g^2d* t])/(Sqrt[4π(σ^2+d t)])))Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;P2[x_,t_]=((ϵ0)/2)Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])-g^2d t/2]*Exp[ϵ0(x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])+ϵ0^2d/2t](-RealSign[(((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0])/(2d t/2))+ϵ0]+Erf[(((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0]+2ϵ0 d/2 t)/(2Sqrt[d t/2]))])/.RealSign[(((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0])/(2d t/2))+ϵ0]->1;
(*,Assumptions->Element[d,Positive]&&Element[g,Positive]&&Element[t,Positive]&&Element[x0,Positive]];*)
b=(( k0-d g)/(2(d/2)));
P3[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(RealAbs[x]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((RealAbs[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((RealAbs[x]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]))/.RealSign[((RealAbs[x]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1;
P3aplus[x_,t_]=P3[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->x};
P3aminus[x_,t_]=P3[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->2a-x};
P3plus[x_,t_]=P3[x,t]/.{RealAbs[x]->x};
P3minus[x_,t_]=P3[x,t]/.{RealAbs[x]->-x};
Print["Solubility check of P3"]
Assuming[x≠ a,Simplify[d/2 D[P3[x,t],{x,2}]+d g* D[P3[x,t]*RealSign[x-a],x]-D[P3[x,t],t]]]
Print["Continuity check of P3 a, 0"]
Simplify[P3aplus[a,t]-P3aminus[a,t]]
Simplify[P3plus[0,t]-P3minus[0,t]]
Print["Boundary check for P3"]
Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&((k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[P3plus[x,t],{x,1}]+d g*P3plus[x,t]-d/2 D[P3minus[x,t],{x,1}]-d g*P3minus[x,t]*(-1)-(k0)P3[a,t]/.x->a]]

Print["Our solution is"]
P[x_,t_]=P1[x,t]+P2[x,t]+P3[x,t]
(*Here we start assuming a to be negative*)
Pplus[x_,t_]=P[x,t]/.{RealAbs[x]->x};
Pminus[x_,t_]=P[x,t]/.{RealAbs[x]->-x};
Paplus[x_,t_]=P[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->x};
Paminus[x_,t_]=P[x,t]/.{x*RealSign[x-a]+2a*UnitStep[-x+a]->2a-x};
Print["Continuity check for the total solution at a, 0"]
Simplify[Paplus[a,t]-Paminus[a,t]]
Simplify[Pplus[0,t]-Pminus[0,t]]
Print["Boundary check at x=a leads to differential equation"]
If[Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&(( k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[Paplus[x,t],{x,1}]+d g*Paplus[x,t]-d/2 D[Paminus[x,t],{x,1}]-d g*Paminus[x,t]*(-1)/.x->a]]!=0,Print["THe problem is you are contraining the boundary conditions one by one, that is why the problem is"],Print["THe method works"]]
Print["Boundary at x=0 leads to differential equation"]
If[Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&(( k0 t+x0)/(d t))≠0&&a<0,Simplify[d/2 D[Pplus[x,t],{x,1}]+d g*Pplus[x,t]-d/2 D[Pminus[x,t],{x,1}]-d g*Pminus[x,t]*(-1)-( k0)Pplus[0,t]/.x->0]]!=0,Print["THe problem is you are contraining the boundary conditions one by one, that is why the problem is"],Print["THe method works"]]
Beep[]

So let us add both the boundary conditions simultaneously, there fore the new PDE could be,

Verifying the result of absolute value potential with the numerical test

Clear[t]
Protect[p1]
ClearAll["Global'*"];
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
x0=5;
σ=1/Sqrt[400];
de=.1;
k=.1;
k0=1;
dThick=40;
dTime=100;
eqSimplee=D[Pe[x,t],t]==de/2*D[Pe[x,t],{x,2}]+k*D[x*Pe[x,t],x]-k0*S[x]*Pe[x,t];
β=1/(4σ^2);
eqInie=Pe[x,0]==√(β/π)*Exp[-β*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
h=-.4;
y=-1.5;
ys=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}]
ab1=Plot[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick}],{t,0,100},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
Show[ab1,p1]

ClearAll["Global`*"];
b=(( k0-d g)/(2(d/2)));
P3[x_,t_]=Exp[-g*(-RealAbs[x0]+x*RealSign[x-a]+2a*UnitStep[-x+a])](((Exp[-(((x+x0)^2)/(4 d t/2))-g^2d t/2])/(Sqrt[4π d t/2]))+(b/2)*Exp[b(x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])+b^2(σ^2+d/2 t)-g^2d t/2]*(-RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]+Erf[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0]+2b(σ^2+d t/2))/(2Sqrt[(σ^2+d t/2)]))]))/.RealSign[((x*RealSign[x-a]+2a*UnitStep[-x+a]+RealAbs[x0])/(2(σ^2+d/2 t)))+b]->1;

ClearAll["Global`*"];
eqn=(d/2)D[P[x,t],{x,2}]+2d g P[x,t]DiracDelta[x-a]+(2d g-k0)DiracDelta[x]P[x,t]==D[P[x,t],t]//TraditionalForm(*eq.1*)
P[x,0]=Exp[g*RealAbs[x0]-g^2d t/2]DiracDelta[x-(Subscript[x, 0])]//TraditionalForm
y1[x_]=1*Exp[I k x]+b*Exp[-I k x];
y2[x_]=c Exp[I k x]+d Exp[-I k x];
y3[x_]=e Exp[I k x];
y1[-a]==y2[-a]
y2[a]==y3[a]
Simplify[y2'[x]-y1'[x]]/.x->-a
Simplify[y3'[a]-y2'[a]]
(*{b1,c1,d1,e1}={b,c,d,e}/.Flatten[Solve[{y1[-a]==y2[-a],y2[a]==y3[a],y2'[-a]-y1'[-a]+2d g y2[-a]==0,y3'[a]-y2'[a]-(k0-2d g)y2[a]==0},{b,c,d,e}]]
{b,c,d,e}={b,c,d,e}/.Flatten[Solve[{y1[0]==y2[0],y2[a]==y3[a],y2'[0]-y1'[0]==-2d g y2[0],y3'[a]-y2'[a]==(k0-2d g)y2[a]},{b,c,d,e}]]
Simplify[b-b1]
Simplify[c-c1]
Simplify[d-d1]
Simplify[e-e1]*)
c=((2I k Exp[-I k a])/(2I k Exp[-I k a]-k1 Exp[-I k a]-((k2 k1*Exp[I 3 k a])/(2I k-k2))));
e=((2I k)/(2I k-k2))c;
d=(e-c)Exp[I 2 k a];
b=(c Exp[-I k a]+d Exp[I k a]-Exp[-I k a])Exp[-I k a];
y1[-a]===y2[-a]
Simplify[y2[a]-y3[a]]
Simplify[(y2'[-a]-y1'[-a])-k1*y1[-a]]
Simplify[y3'[a]-y2'[a]-k2*y2[a]]
Simplify[b]
Simplify[c]
Simplify[d]
Simplify[e]

Assuming[((d g t-RealAbs[x0]-2 a UnitStep[a])/(d t))≠0&&a≠0&&x0>0&&d∈Positive&&g∈Positive&&(( k0 t+x0)/(d t))≠0&&a<0,Simplify[%2174]]

ClearAll["Global`*"]
a=0;
y1[x_]=a1*Exp[I (p+q) x]+b1*Exp[I (-p+q) x];
y2[x_]=c1*Exp[I (p-q) x]+d1 Exp[I (-p-q) x];
a1=1;
d1=1;
{{b1,c1}}={b1,c1}/.Simplify[Solve[{y1[a]==y2[a],d(y2'[a]-y1'[a])==-2d g*y2[0]},{b1,c1}]]
Simplify[b1]
Simplify[c1]
q=I g/2;
p=q+Sqrt[q^2+4s/d]

Together[p]

b

ClearAll["Global`*"]
y1[x_]=a*Exp[I k x]+b Exp[-I k x];
y2[x_]=c Exp[I k x];
p0[x_]=Exp[I k x]-0(((E^(2 I a k)) g)/(g+I k))Exp[I k Abs[x]];
a=1;
{{b,c}}={b,c}/.Simplify[Solve[{y1[0]==y2[0],(d/2)(y2'[0]-y1'[0])==k0(p0[0]+0y2[0])},{b,c}]]
Clear[b,c];
{{b,c}}={b,c}/.Simplify[Solve[{y1[0]==y2[0],(d/2)(y2'[0]-y1'[0])==k0(y2[0])},{b,c}]]

ClearAll["Global`*"]
y1[x_]=a1*Exp[I k x]+b Exp[-I k x];
y2[x_]=c Exp[I k x];
p0[x_]=Exp[I k x]+0(((-I g)/(k+I g)))Exp[I k Abs[x]];
a1=1;
{{b,c}}={b,c}/.Simplify[Solve[{y1[a]==y2[a],(d/2)(y2'[a]-y1'[a])==-2d g*y2[a]+d g*y2[a]},{b,c}]]

An attempt to do the procedure using mathematica;

Integrate[((ϵ0)/2)*Exp[ϵ0(Abs[x]+x0)+ϵ0^2(σ^2+d/(4k)*Exp[2 k t])+k t]*(-Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]+Erf[((Abs[x]+x0+ϵ0(σ^2+d/(4k)*Exp[2 k t]))/(2Sqrt[(σ^2+d/(4k)*Exp[2 k t])]))]),{x,-∞,x1*Exp[k t]}]

ClearAll[x]
F[x_,t_]:=Simplify[(1/2) (Sqrt[π]) ((1/(Sqrt[(k/(d (E^(2 k t))+4 k (σ^2)))]))+((((E^(k t)) x+x0) Erf[(Sqrt[((k ((((E^(k t)) x+x0))^2))/(d (E^(2 k t))+4 k (σ^2)))])])/(Sqrt[((k ((((E^(k t)) x+x0))^2))/(d (E^(2 k t))+4 k (σ^2)))])))((Exp[k t])/(2Sqrt[π]Sqrt[σ^2+(d/(4k))Exp[2k t]]))+ϵ0*Exp[ϵ0(σ^2+(d/(4k))*Exp[2k t])]((Exp[ϵ0 x Exp[k t]])/(ϵ0))(-(*Sign[((x*Exp[k t]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]*)1+Erf[(x*Exp[k t]+x0+ϵ0 2(σ^2+d/(4k)*Exp[2 k t]))/(2Sqrt[(σ^2+d/(4k)*Exp[2 k t])])])+ϵ0*Exp[-ϵ0^2(σ^2+(d/(4k))Exp[2k t])]((2/(ϵ0)))*UnitStep[-x0-2ϵ0(σ^2+(d/(4k))Exp[2k t])]+ϵ0*Exp[ϵ0^2(σ^2+(d/(4k))Exp[2k t])+ϵ0 x0](((-1)/(ϵ0)))*UnitStep[-x0-2ϵ0(σ^2+(d/(4k))Exp[2k t])]*Exp[((-(x0+2ϵ0(σ^2+(d/(4k))Exp[2k t]))^2)/(4(σ^2+(d/(4k))Exp[2k t])))]*Exp[((x0^2)/(4(σ^2+(d/(4k))Exp[2k t])))](Erf[((x Exp[k t]+x0)/(2Sqrt[(σ^2+(d/(4k))Exp[2k t])]))]-Erf[(x0/(2Sqrt[(σ^2+(d/(4k))Exp[2k t])]))])]
P[x_,t_]=(∂_x)F[x,t]
ϵ0=k0/d;
a=d/2*D[P[x,t],{x,1}]+k x P[x,t]/.x->0
b=-d/2*D[P[x,t],{x,1}]+k x P[x,t]/.x->0
c=-k0 Exp[k t]*P[0,t]
Simplify[a+b+c]

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]

Let us cal the alternative method for help;

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
P[x_,t_]=Simplify[((Exp[((-(x+x0)^2)/(4(σ^2+d*Exp[2 k t]/(4k))))+k t])/(Sqrt[(σ^2+d*Exp[2 k t]/(4k))4π]))+((ϵ0)/2)*Exp[k t+ϵ0(Abs[x]+x0)+ϵ0^2(σ^2+d/(4k)*Exp[2 k t])]*(-(*Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]*)1+Erf[((Abs[x]+x0+2ϵ0(σ^2+d/(4k)*Exp[2 k t]))/(2Sqrt[(σ^2+d/(4k)*Exp[2 k t])]))])/.{x->x*Exp[k t],d*Exp[2 k t]->d*(-1+Exp[2k t])},Element[k,Reals]&&Element[t,Reals]]
ϵ0=(k0/d);
Peplus[x_,t_]=P[x,t]/.Abs[x]->x;
Peminus[x_,t_]=P[x,t]/.Abs[x]->-x;
Print["Analytical check"]
Simplify[d/2*D[Peplus[x,t],{x,2}]+k*D[x*Peplus[x,t],x]-D[Peplus[x,t],t],Assumptions->Element["Global`*",Reals]]
Simplify[d/2*D[Peplus[x,t],{x,2}]+k*D[x*Peplus[x,t],x]-D[Peplus[x,t],t],Assumptions->Element["Global`*",Reals]]
a=d/2*D[Peplus[x,t],{x,1}]+k x Peplus[x,t]/.x->0;
b=-d/2*D[Peminus[x,t],{x,1}]-k x Peminus[x,t]/.x->0;
c=-k0 Exp[k t] P[0,t];
Print["Boundary condition check"]
Simplify[a+b+c]

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
P[x_,t_]=Simplify[((Exp[((-(x+x0)^2)/(4(σ^2+d*Exp[2 k t]/(4k))))-2k t])/(Sqrt[(σ^2+d*Exp[2 k t]/(4k))4π]))+((ϵ0)/2)*Exp[-2k t+ϵ0(Abs[x]+x0)+ϵ0^2(σ^2+d/(4k)*Exp[2 k t])]*(-(*Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]*)1+Erf[((Abs[x]+x0+2ϵ0(σ^2+d/(4k)*Exp[2 k t]))/(2Sqrt[(σ^2+d/(4k)*Exp[2 k t])]))])/.{x0->x0*Exp[-k t],d*Exp[2 k t]->d*(1-Exp[-2k t])},Element[k,Reals]&&Element[t,Reals]]
ϵ0=(k0/d);
σ=0;
Peplus[x_,t_]=P[x,t]/.Abs[x]->x;
Peminus[x_,t_]=P[x,t]/.Abs[x]->-x;
Print["Analytical check"]
Simplify[d/2*D[Peplus[x,t],{x,2}]+k*D[x*Peplus[x,t],x]-D[Peplus[x,t],t],Assumptions->Element["Global`*",Reals]]
Simplify[d/2*D[Peplus[x,t],{x,2}]+k*D[x*Peplus[x,t],x]-D[Peplus[x,t],t],Assumptions->Element["Global`*",Reals]]
a=d/2*D[Peplus[x,t],{x,1}]+k x Peplus[x,t]/.x->0;
b=-d/2*D[Peminus[x,t],{x,1}]-k x Peminus[x,t]/.x->0;
c=-k0 Exp[k t] P[0,t];
Print["Boundary condition check"]
Simplify[a+b+c]

ClearAll["Global`*"]
Unprotect["Global`*"]
Remove["Global`*"]
P[x_,t_]=Simplify[((Exp[((-(x+x0)^2)/(4(σ^2+d*Exp[2 k t]/(4k))))+k t])/(Sqrt[(σ^2+d*Exp[2 k t]/(4k))4π]))+((ϵ0)/2)*Exp[ϵ0(Abs[x]+x0)+ϵ0^2(σ^2+d/(4k)*Exp[2 k t])+k t]*(-1+Erf[((Abs[x]+x0+2ϵ0(σ^2+d/(4k)*Exp[2 k t]))/(2Sqrt[(σ^2+d/(4k)*Exp[2 k t])]))])/.{x->x*Exp[k t],d*Exp[2 k t]->d*(-1+Exp[2k t])},Element[k,Reals]&&Element[t,Reals]];
(*Sign[((Abs[x]+x0)/(2(σ^2+d/(4k)*Exp[2 k t])))+ϵ0]*)
P[x,t]
k0=10;
ϵ0=(k0/d);
a=ϵ0;
σ=0;
P[x,t]
k=1;
d=.3;
x0=5;
p1=Plot[{P[x,.1],P[x,.3],P[x,1],P[x,2],P[x,5],P[x,10]},{x,-20,20},PlotRange->Full]
Print["When k0 is 10"]
s1=P[x,t]
ClearAll[k0,ϵ0,a]
k0=0;
ϵ0=(k0/d);
a=ϵ0;
Print["When k0 is 0"]
s2=P[x,t]
p2=Plot[{P[x,.1],P[x,.3],P[x,1],P[x,2],P[x,5],P[x,10]},{x,-20,20},PlotRange->Full]
u[x_,t_]=Sqrt[(k/(2π*d/2(1-Exp[-2k t])))]*Exp[((-k)/(2d/2))*(((x+x0*Exp[-k t])^2)/((1-Exp[-2k t])))];
p3=Plot[{u[x,.1],u[x,.3],u[x,1],u[x,2],u[x,5],u[x,10]},{x,-20,20},PlotRange->Full]
Show[p1,p2,p3]
Print["Difference beween k0=10 and k0=0 case"]
Simplify[s1-s2]
Print["Difference beween k0=0 and u[x,t]"]
Simplify[s2-u[x,t]]
Print["Difference beween k0=10 and u[x,t]"]
Simplify[s1-u[x,t]]

(*
t={0.1,5,20};
NIntegrate[{P[x,t]},{x,-20,20}]*)
(*Q[t_]=1+Exp[k t+ϵ0*x0+ϵ0^2(σ^2+d*Exp[2 k t]/(4k))](Sign[(x0/((2(σ^2+d*Exp[2 k t]/(4k)))))+a]-Erf[((x0+2a(σ^2+ d*Exp[2 k t]/(4k)))/((2Sqrt[(σ^2+ d*Exp[2 k t]/(4k))])))])
+2*UnitStep[-x0-2a(σ^2+ d*Exp[2 k t]/(4k)t)]*Exp[-ϵ0^2(σ^2+ d*Exp[2 k t]/(4k) t)]-Erfc[(x0/((Sqrt[4(σ^2+d*Exp[2 k t]/(4k) t)])))]/.{x->x*Exp[k t],d*Exp[2 k t]->d*(-1+Exp[2k t])};
Plot[{Q[t],k0 Exp[k t] P[0,t]},{t,0,100},PlotRange->Full]*)

Plot3D[%501,{t,0,100},{x,-30,30},PlotRange->All]

ClearAll["Global`*"];
Integrate[Exp[((-(x*Exp[k t]+x0)^2)/(4(σ^2+(d/(4k))(Exp[2 k t]-1))))],{x,0,∞}]

Clear[t]
Protect[p1]
ClearAll["Global'*"];
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
σ=1/Sqrt[40];
k0=3;
k=1;
d=.5;
x0=3;
dThick=40;
dTime=100;
eqSimplee=D[Pe[x,t],t]==d/2*D[Pe[x,t],{x,2}]+k*D[x*Pe[x,t],x]-k0*Exp[k t]*S[x]*Pe[x,t];
β=1/(4σ^2);
eqInie=Pe[x,0]==√(β/π)*Exp[-β*(x+x0)^2];
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pe[dThick,t]==0; 
h=-.4;
y=-1.5;
ys=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}]
Plot3D[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick},{t,0,50},PlotRange->Full]
Table[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick}],{t,0,20}]
ListPlot[%]
(*ab1=Plot[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick}],{t,0,100,5},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
Show[ab1,p1]*)

Flatten[%717]

ListPlot[%719]

The kenkre’s method of finding the answer to arbitrarily placed sink.

ClearAll["Global`*"];
p1[x_]=a1*Exp[(1+Sqrt[1+((4l ϵ)/Γ)])(x/2l)];(*x<x0*)
p2[x_]=b*Exp[(1+Sqrt[1+((4l ϵ)/Γ)])(x/2l)]+c*Exp[(1-Sqrt[1+((4l ϵ)/Γ)])(x/2l)];(*x>x0&&x<a*)
p3[x_]=d*Exp[(-1-Sqrt[1+((4l ϵ)/Γ)])(x/2l)];(*x>a*)
eqn1=p1[x0]==p2[x0];
eqn2=p2'[x0]-p1'[x0]-p2[x0]+p1[x0]==k0*p2[x0];
eqn3=p2[a]==p3[a];
eqn4=p3'[a]-p2'[a]+p3[a]+p2[a]==0;
{a1,b,c,d}={a1,b,c,d}/.Flatten[Solve[{p1[x0]==p2[x0],d(p2'[x0]-p1'[x0])+Γ(-p2[x0]+p1[x0])==-1,p2[a]==p3[a],d(p3'[a]-p2'[a])+Γ(p3[a]+p2[a])==0},{a1,b,c,d}]]
Clear[a];
a
p1[x]
Print["With a=0"]
a=0;
p1[x]
Clear[a]
a
Simplify[p2[x]]
Print["With a=0"]
a=0;
Simplify[p2[x]]
Clear[a]
a
Simplify[p3[x]]
Print["With a=0"]
a=0;
Simplify[p3[x]]

Taking Inverse Laplace transform of the kenkre solution

ClearAll["Global`*"];
P[x_,x0_,s_]=((Exp[-(((Abs[x]-Abs[x0]))/(2l))])/(Γ Sqrt[1+((4 l s)/Γ)]))(Exp[-Sqrt[1+((4 l s)/Γ)]((Abs[x-x0])/(2l))]+0((Exp[-Sqrt[1+((4l s)/Γ)](((Abs[x]+Abs[x0]))/(2l))])/(Sqrt[1+((4 l s)/Γ)]-1)));
InverseLaplaceTransform[P[x,x0,s],s,t]

THe initial condition test
ClearAll["Global`*"];P[x_,t_]=(Exp[-g*(x*RealSign[x-a]+2a*UnitStep[-x+a]-x0*RealSign[x0+a]-2a UnitStep[x0+a])-g^2d* t]/Sqrt[4π(σ^2+d t)])Exp[-(((x+x0)^2)/(4(σ^2+d t)))]/.t->t/2;
P[x,0]
x=-x0;
Simplify[P[x,0]]
Clear[x]
Assuming[x>a&&x0>-a,Simplify[(x*RealSign[x-a]+2a*UnitStep[-x+a]-x0*RealSign[x0+a]-2a UnitStep[x0+a])]]
Assuming[x<a&&x0>-a,Simplify[(x*RealSign[x-a]+2a*UnitStep[-x+a]-x0*RealSign[x0+a]-2a UnitStep[x0+a])]]
Assuming[x>a&&x0<-a,Simplify[(x*RealSign[x-a]+2a*UnitStep[-x+a]-x0*RealSign[x0+a]-2a UnitStep[x0+a])]]
Assuming[x<a&&x0<-a,Simplify[(x*RealSign[x-a]+2a*UnitStep[-x+a]-x0*RealSign[x0+a]-2a UnitStep[x0+a])]]
((E^(-((((x+x0))^2)/(4 (σ^2)))-g (x RealSign[-a+x]-x0 RealSign[a+x0]+2 a UnitStep[a-x]-2 a UnitStep[a+x0])))/(2 (Sqrt[π]) (Sqrt[(σ^2)])))
(1/(2 (Sqrt[π]) (Sqrt[(σ^2)])))
-2 a+x-x0
-x-x0
x+x0
2 a-x+x0
\exp\left(-g (-\text{x0}\text{sgn}(-a-\text{x0})+2 a\theta (a+\text{x0})+\left|
\text{x0}\right|)-d g^2 t\right)
TeXForm[(Exp[-g*((x*RealSign[x-a]+2a*UnitStep[-x+a])+RealAbs[x0])-g^2d*t]/Sqrt[4π(σ^2+d t)])Exp[-((x+x0)^2/(4(σ^2+d t)))]]
\frac{\exp\left(-g (x\text{sgn}(x-a)+2 a\theta (a-x)+\left|
\text{x0}\right|)-d g^2 t-\frac{(x+\text{x0})^2}{4\left(d t+\sigma^2\right)}\right)}{2\sqrt{\pi}\sqrt{d t+\sigma^2}}

ClearAll["Global`*"];Pe[x_,t_]=((1/(2Sqrt[π(σ^2+ d t)])))Exp[-(((x+x0)^2)/(4 (σ^2+ d t)))]+(k2/(4d))Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x*Sign[x-b]+2b*Unit[b-x]+x0)](-RealSign[((x*Sign[x-b]+2b*Unit[b-x]+x0)/((2(σ^2+ d t))))+2ϵ2]+Erf[((x*Sign[x-b]+2b*Unit[b-x]+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+Heaviside[a-x]((k1/(4d))Exp[ϵ2^2*(σ^2+ d t)+ϵ2(-x+2b+x0)](-RealSign[((-x+2b+x0)/((2(σ^2+ d t))))+2ϵ2]+Erf[((-x+2b+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])-(k1/(4d))Exp[ϵ2^2*(σ^2+ d t)+ϵ2(x+2(b-a)+x0)](-RealSign[((x+2(b-a)+x0)/((2(σ^2+ d t))))+2ϵ2]+Erf[((x+2(b-a)+x0+2ϵ2(σ^2+ d t))/(2Sqrt[(σ^2+ d t)]))])+((ϵ1)/2)(Erfc[((x+x0)/(2Sqrt[(σ^2+ d t)]))]+Erfc[((-x+2b+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((x+2(b-a)+x0)/(2Sqrt[(σ^2+ d t)]))]-Erfc[((-x+2a+x0)/(2Sqrt[(σ^2+ d t)]))]));
TeXForm[Pe[x,t]]