The paper presents all the results of “Exact time-domain solution of the SchrÖdinger equation for a new scattering model”

The paper presents a new method for solving a time-dependent Schrodinger equation with an ultrashort barrier. So far, there are no barriers other than the Dirac delta barrier has an solution. THis is   a first time a barrier with a short range has an exact solution

The following are the some of the useful integrals;

ClearAll["Global`*"];
Integrate[Exp[-a u^2]/(u+c+I b),{u,-∞,∞},Assumptions->b<0&&c>0&&a∈Positive]
Integrate[Exp[-a u^2]/(u+c+I b),{u,-∞,∞},Assumptions->b>0&&c>0&&a∈Positive]
Integrate[Exp[-a u^2]/(u+c+I b),{u,-∞,∞},Assumptions->b>0&&c<0&&a∈Positive]
Integrate[Exp[-a u^2]/(u+c+I b),{u,-∞,∞},Assumptions->b<0&&c<0&&a∈Positive]

THis section presents the result of a Dirac delta barrier;

Let us visualize the result of Fei sun;
THe derived expression matches with that of the Fei sun. A small variation in the magnitudes can be due to numerical inaccuracy of his code

Unprotect["Global`*"];
Remove["Global`*"];
ClearAll["Global`*"]
β0=((m V0)/(ℏ^2));
σ=Sqrt[2];
V0=10;
ℏ=1.055;
m=9.109;
k0=3;
x0=8;
ψimage[x_,t_]:=(((σ^2)/(2 π (σ^2+I ℏ t/(2m)))))^.25 Exp[((-(x+x0)^2)/(4 (σ^2+I ℏ t/(2m))))+((I σ^2(k0(x+x0)-((ℏ k0^2)/(2m ))t))/((σ^2+I ℏ t/(2m))))]-(((σ^2)/(2 π (σ^2+I ℏ t/(2m)))))^.25 Exp[((-(RealAbs[x]+x0)^2)/(4 (σ^2+I ℏ t/(2m))))+((I σ^2(k0(RealAbs[x]+x0)-((ℏ k0^2)/(2m ))t))/((σ^2+I ℏ t/(2m))))];
ψ[x_,t_]:=(((σ^2)/(2 π (σ^2+I ℏ t/(2m)))))^.25 Exp[((-(x+x0)^2)/(4 (σ^2+I ℏ t/(2m))))+((I σ^2(k0(x+x0)-((ℏ k0^2)/(2m ))t))/((σ^2+I ℏ t/(2m))))]-((I m V0)/(ℏ^2))Sqrt[(((σ^2+I ℏ t/(2m))/π))](((σ^2)/(2 π (σ^2+I ℏ t/(2m)))))^.25 Exp[((-(RealAbs[x]+x0)^2)/(4 (σ^2+I ℏ t/(2m))))+((I σ^2(k0(RealAbs[x]+x0)-((ℏ k0^2)/(2m ))t))/((σ^2+I ℏ t/(2m))))](-I π Exp[(σ^2+I ℏ t/(2m))(- I κ+β0)^2])(Erfc[Sqrt[(σ^2+I ℏ t/(2m))](-I κ+ β0)]);
κ=((k0 σ^2)/((σ^2+I ℏ t/(2m))))+((I (RealAbs[x]+x0))/(2(σ^2+I ℏ t/(2m))));
Plot3D[Abs[ψ[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψ[x,t]]^2,{x,-30,30},{t,0,50}]
t=25;
p2=Plot[{Abs[ψ[x,t]]^2},{x,-8,8},PlotTheme->"Scientific",PlotRange->Full,PlotLegends->Placed[{Style["DDP",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}],PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]}]
p5=Plot[{-Abs[ψ[x,t]]^2+Abs[ψimage[x,t]]^2},{x,-8,8},PlotLegends->Automatic,PlotRange->Full]
Clear[t];
t=10;
p1=Plot[{Abs[ψ[x,t]]^2,Abs[ψimage[x,t]]^2},{x,-100,100},PlotLegends->Automatic,PlotRange->Full];
Clear[t];
t=110;
p3=Plot[{Abs[ψ[x,t]]^2,Abs[ψimage[x,t]]^2},{x,-100,100},PlotLegends->Automatic,PlotRange->Full];
NIntegrate[Abs[ψ[x,1]]^2,{x,-100,100}](*sink free N=1.0571880110215248`*)
NIntegrate[Abs[ψ[x,10]]^2,{x,-100,100}](*sink free N=2.2581008643542058`*)
NIntegrate[Abs[ψ[x,100]]^2,{x,-1000,1000}](*sink free N=7.071774812605344`*)

The numerical verification for the analytical solution

Numerical evaluation of the wave-packet using ND solve of Mathematica

Unprotect["Global`*"];
ClearAll["Global`*"];
Remove["Global`*"];NDSolve[{I D[u[t,x],t]+D[u[t,x],x,x]==1/(1+Exp[-x])*u[t,x],u[0,x]==DiracDelta[x+5],u[t,-100]==u[t,100]},u,{t,0,10},{x,-100,100}]
t={1,4,10};
aa1=Plot[Abs[%%[[1,1,2]][t,x]]^2,{x,-20,20},PlotPoints->100]
ab1=Plot[Abs[((1/(Sqrt[2]π (I t))))^.5 Exp[((-(x+5)^2)/(2 (I t)))]]^2,{x,-20,20},PlotRange->Full]
Show[ab1,aa1]

Unprotect["Global'*"];
ClearAll["Global'*"];
Remove["Global`*"];
x0=1;
σ=1/Sqrt[40];
m=1;
ℏ=1;
V0=1;
d=-.7;
dThick=40;
dTime=100;
eqSimplee=I*ℏ*D[Psi[x,t],t]==-((ℏ^2)/(2m))*D[Psi[x,t],{x,2}]+0*V0*UnitStep[x+d]*UnitStep[d-x]*Psi[x,t];
β=1/(4σ^2);
eqInie=Psi[x,0]==((1/(2π σ^2)))^0.25*Exp[-β*(x+x0)^2];
eqBlefte=Abs[Psi[-dThick,t]]==0;
eqBrighte=Abs[Psi[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,eqBrighte},Psi[x,t],{x,-dThick,dThick},{t,0,dTime}]
Plot3D[Re[Psi[x,t]/.res],{x,-dThick,dThick},{t,0,dTime},PlotRange->Full]
t=10;
aa1=Plot[Evaluate[Abs[Psi[x,t]/.res]^2],{x,-dThick,dThick},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, g])(t)",14]},PlotStyle->{Dashed}]
ab1=Plot[Abs[(((σ^2)/(2 π (σ^2+I ℏ t/(2m)))))^.25 Exp[((-(x+x0)^2)/(4 (σ^2+I ℏ t/(2m))))]]^2,{x,-dThick,dThick}]
Show[aa1,ab1]

This section presents the space-time dynamics of the wave packet interacting with the barrier. The effect of barrier height is shown and the result is compared with that of the delta barrier

Let us introduce the wavefunction of the Gaussian intial wave packet as ψg; V0 effect graphs are here; t effect graphs next cell

Protect[p2];
ClearAll["Global`*"]
Remove["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=Sqrt[2];
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
d=.5;
k0=3;
ν0=0;
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(2k0*(σ^2)+I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
ψg[x_,t_]=HeavisideTheta[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi-UnitStep[x-d]2d)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))-((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(xi-d)-((ℏ k0^2t)/(2m)))-(((-d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+(((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2-I*ν0/ℏ]*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[1])])+I*π*Exp[(-I f[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[-1])])))];
V0=1;
Print["The pole value is"]
k[1]
Plot3D[Abs[ψg[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψg[x,t]]^2,{x,-30,30},{t,0,50}]
Plot[(1/(1^2+(((m*ℏ (En-V0)2d)/(ℏ^2 *Sqrt[2 m En])))^2)),{En,0,30},PlotLegends->"T(E)"]
t=25;
p1=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotStyle->{RGBColor[0,0,1],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["(Subscript[V, 0])=1",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["(Subscript[V, 0])=1",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}]]
Clear[V0]
V0=10
p2a=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotStyle->{RGBColor[0,0,1],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["(Subscript[V, 0])=10",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["(Subscript[V, 0])=10",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}]]
Clear[V0]
V0=100;
p3=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotStyle->{RGBColor[0,0,0],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["(Subscript[V, 0])=100",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["(Subscript[V, 0])=100",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}]]
plots=Show[p1,p2,p3]
Style[plots,ImageSize->1]
Clear[V0,t,k];
Remove[k[i]]
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
V0=100;
t=25;
Print["The pole value is"]
k[1]
p4=Plot[{Abs[ψg[x,t]]^2},{x,-10,20},PlotRange->Full,PlotStyle->{RGBColor[0,0,1],Dotted}]
Show[p2,p2a]
Clear[t]
NIntegrate[Abs[ψg[x,1]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,10]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,50]]^2,{x,-100,100},MaxRecursion->20]

The dynamics of the probability distribution is calculated.

Protect[p2];
ClearAll["Global`*"]
Remove["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=Sqrt[2];
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
d=.5;
k0=3;
ν0=0;
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(2k0*(σ^2)+I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
ψg[x_,t_]=HeavisideTheta[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi-UnitStep[x-d]2d)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))-((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(xi-d)-((ℏ k0^2t)/(2m)))-(((-d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+(((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2-I*ν0/ℏ]*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[1])])+I*π*Exp[(-I f[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[-1])])))];
V0=5;
Print["The pole value is"]
k[1]
Plot3D[Abs[ψg[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψg[x,t]]^2,{x,-30,30},{t,0,50}]
Plot[(1/(1^2+(((m*ℏ (En-V0)2d)/(ℏ^2 *Sqrt[2 m En])))^2)),{En,0,30},PlotLegends->"T(E)"]
t=5;
p1=Plot[{Abs[ψg[x,t]]^2},{x,-10,8},PlotStyle->{RGBColor[0,0,1],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=5",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+1}}]]
Clear[t]
t=10
p2a=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotStyle->{RGBColor[0,0,0],Dotted},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=10",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}]]
Clear[t]
t=25;
p3=Plot[{Abs[ψg[x,t]]^2},{x,-20,20},PlotStyle->{RGBColor[0,0,1],Dashed},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=25",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[t]
t=50;
p4=Plot[{Abs[ψg[x,t]]^2},{x,-20,20},PlotStyle->{RGBColor[1,1,0],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=50",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5-.5}}]]
plots=Show[p3,p2a,p4,p1]
Style[plots,ImageSize->1]
NIntegrate[Abs[ψg[x,1]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,10]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,50]]^2,{x,-100,100},MaxRecursion->20]

In this section, we calculate various properties of the wavepacket and the barrier. First we calculated the time-dependent dispersion of the wavepacket. Next, we calculate the trajectory of the  mean-position of the wavepacket.  Next the dwell probability inside the barrier, reflection probability and transmission probability are calculated together and is presented.

Computing the probability properties and the statistical properties of the wavepacket

ClearAll["Global`*"]
Remove["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=Sqrt[2];
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
k0=3;
ν0=0;
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(2k0*(σ^2)+I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
ψg[x_,t_]=HeavisideTheta[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi-UnitStep[x-d]2d)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))-((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(xi-d)-((ℏ k0^2t)/(2m)))-(((-d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+(((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2-I*ν0/ℏ]*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[1])])+I*π*Exp[(-I f[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[-1])])))];
d=.5;
V0=1;
var1=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p0=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["a",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]];
Clear[σ]
σ=Sqrt[3];
var1a=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p0a=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["σ=√3",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]];
Clear[k0];
k0=4;
var1b=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p0b=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["(Subscript[k, 0])=4",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]];
Clear[σ]
σ=2;
var1c=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p0c=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["σ=2",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]];
Clear[d,V0,σ,k0]
d=.5;
V0=2;
k0=3;
σ=Sqrt[2];
var2=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p01=ListLinePlot[%,PlotStyle->{Dotted},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["b",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-1}}]];
Clear[d,V0]
d=.5;
V0=3;
var3=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p02=ListLinePlot[%,PlotStyle->{Dashed},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["c",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-2}}]];
Clear[d,V0]
d=.2;
V0=3;
var4=Table[{t,Abs[Sqrt[(NIntegrate[x ^2Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]-(NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20])^2)/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]]]},{t,1,50,2}];
p03=ListLinePlot[%,PlotStyle->{Gray},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|(Subscript[σ, t])|",14]},PlotLegends->Placed[{Style["d",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-3}}]];
Show[p0,p01,p02,p03,p0a,p0c,p0b]

Abs[Sqrt[((sec1[[All,2]]-mean1[[All,2]]^2)/(norm1[[All,2]]))]]

ListPlot[%486]

ClearAll["Global`*"]
Remove["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=Sqrt[2];
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
k0=3;
ν0=0;
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(2k0*(σ^2)+I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
ψg[x_,t_]=HeavisideTheta[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi-UnitStep[x-d]2d)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))-((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(xi-d)-((ℏ k0^2t)/(2m)))-(((-d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+(((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2-I*ν0/ℏ]*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[1])])+I*π*Exp[(-I f[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[-1])])))];
d=.5;
V0=1;
xmean1=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p0=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["a",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[σ]
σ=10;
xmean1a=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p0a=ListLinePlot[%,PlotStyle->{Dotted,Blue},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["e",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[k0];
k0=1;
xmean1b=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p0b=ListLinePlot[%,PlotStyle->{Dashed,Blue},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["f",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[σ]
σ=2;
xmean1c=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p0c=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["g",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[d,V0,σ,k0]
d=.5;
V0=2;
k0=3;
σ=Sqrt[2];
xmean2=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p01=ListLinePlot[%,PlotStyle->{Dotted},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["b",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-1}}]]
Clear[d,V0]
d=.5;
V0=3;
xmean3=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p02=ListLinePlot[%,PlotStyle->{Dashed},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["c",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-2}}]]
Clear[d,V0]
d=.2;
V0=3;
xmean4=Table[{t,NIntegrate[x Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]/NIntegrate[Abs[ψg[x,t]]^2,{x,-100,100},MaxRecursion->20]},{t,1,50,2}];
p03=ListLinePlot[%,PlotStyle->{Gray},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["<x(Subscript[>, t])",14]},PlotLegends->Placed[{Style["d",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5-3}}]]
Show[p0,p01,p02,p03,p0a,p0c,p0b]

Show[p01,p03,p02,p0]

Clear[d,V0,n,pplus,p1,pdwell,p2,pminus,p3]
d=.5;
V0=1;
n=NIntegrate[Abs[ψg[x,100]]^2,{x,d,1000},MaxRecursion->20]
pplus=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,d,100},MaxRecursion->20]},{t,1,50,2}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,d,1000},MaxRecursion->20]}],
{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,d,100000},MaxRecursion->20]}]
p1=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(>)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
pdwell=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,-d,d},MaxRecursion->20]},{t,1,50,5}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,-d,d},MaxRecursion->20]}],{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,-d,d},MaxRecursion->20]}]
p2=ListPlot[%,PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(<>)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
pminus=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,-100,-d},MaxRecursion->20]},{t,1,50,2}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,-1000,-d},MaxRecursion->20]}],{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,-100000,-d},MaxRecursion->20]}]
p3=ListLinePlot[%,PlotStyle->{Dashed},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(<)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Show[p1,p2,p3]

This section presents miscellaneous trials;  and all these results are not succesful. Using the correlation function we have tried to match the reflection and transmission coefficient.
the plots for the dynamics of an exponential wavepacket and delta function wavepacket is attempted.

momentum distribution; correlation function in time and energy domain

ClearAll["Global`*"]
Remove["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=Sqrt[2];
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
ν0=0;
k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(2k0*(σ^2)+I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
ψi[x_]=((1/(2 π σ^2)))^0.25 Exp[-(((x+xi)^2)/(4 σ^2))-I k0(x+xi)-I ν0/ℏ];
ψg[x_,t_]=HeavisideTheta[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi-UnitStep[x-d]2d)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))-((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(xi-d)-((ℏ k0^2t)/(2m)))-(((-d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+(((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2-I*ν0/ℏ]*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[1])])+I*π*Exp[(-I f[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f[x,t]+k[-1])])))];
d=.5;
T=30;
V0=5;
dt=1;
wavecorr=Append[Append[Table[{t,Abs[NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(t+T)/ℏ],{x,-100,100},MaxRecursion->20]]},{t,0,10,dt}],{100,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(100+T)/ℏ],{x,-1000,1000},MaxRecursion->20]}],
{1000,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(1000+T)/ℏ],{x,-1000,1000},MaxRecursion->20]}]
p0=ListLinePlot[Abs[%],PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|C(t)|",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
totE=Fourier[Abs[wavecorr[[All,2]]]]
p1=ListLinePlot[Abs[%],PlotStyle->{Thick},PlotTheme->"Scientific",PlotRange->{0,1},FrameLabel->{Style["E",14],Style["σ(E)",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
transcorr=Append[Append[Table[{t,Abs[NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(t+T)/ℏ]*UnitStep[x],{x,-100,100},MaxRecursion->20]]},{t,0,10,dt}],{100,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(100+T)/ℏ]*UnitStep[x],{x,-1000,1000},MaxRecursion->20]}],
{1000,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(1000+T)/ℏ]*UnitStep[x],{x,-1000,1000},MaxRecursion->20]}];
p2=ListLinePlot[Abs[%],PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|C(t)|",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
transE=Abs[Fourier[transcorr[[All,2]]]]/Abs[Fourier[wavecorr[[All,2]]]]
p3=ListLinePlot[%,PlotStyle->{Thick},PlotRange->{0,1},PlotTheme->"Scientific",FrameLabel->{Style["E",14],Style["T(E)",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
refcorr=Append[Append[Table[{t,Abs[NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(t+T)/ℏ]*UnitStep[-x],{x,-100,100},MaxRecursion->20]]},{t,0,10,dt}],{100,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(100+T)/ℏ]*UnitStep[-x],{x,-1000,1000},MaxRecursion->20]}],
{1000,NIntegrate[Conjugate[ψg[x,T/ℏ]]*ψg[x,(1000+T)/ℏ]*UnitStep[-x],{x,-1000,1000},MaxRecursion->20]}];
p4=ListLinePlot[Abs[%],PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|C(t)|",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
refE=Abs[Fourier[refcorr[[All,2]]]]/Abs[Fourier[wavecorr[[All,2]]]]
p5=ListLinePlot[%,PlotStyle->{Thick},PlotRange->{0,1},PlotTheme->"Scientific",FrameLabel->{Style["E",14],Style["T(E)",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]

NIntegrate[Conjugate[ψg[x,0/ℏ]]*ψg[x,0/ℏ],{x,-100,100}]
NIntegrate[Abs[ψg[x,0/ℏ]]^2,{x,-100,100}]
Abs[NIntegrate[Conjugate[ψi[x]]*ψi[x],{x,-100,100}]]
NIntegrate[Abs[ψi[x]]^2,{x,-100,100}]

Abs[0.9420899655084068` -0.5147808275719644` I]

Table[{t,Abs[NIntegrate[Conjugate[ψg[x,30/ℏ]]*ψg[x,(t+30)/ℏ],{x,-100,100},MaxRecursion->20]]},{t,0,2}]
NIntegrate[Abs[ψg[x,30/ℏ]]^2,{x,-100,100}]

Table[{t,Abs[NIntegrate[Conjugate[ψg[x,30/ℏ]]*ψg[x,(t+30)/ℏ],{x,-100,100},MaxRecursion->20]]},{t,0,10}]

p0

wavecorr=Append[Append[Table[{t,NIntegrate[ψi[x]*ψg[x,t],{x,-100,100}]},{t,0,50,5}],{100,NIntegrate[ψi[x]*ψg[x,100],{x,-100,100}]}],
{1000,NIntegrate[ψi[x]*ψg[x,1000],{x,-100,100}]}];
p0=ListLinePlot[Abs[%],PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["|C(t)|",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
totE=Abs[Fourier[wavecorr[[All,2]]]]
p1=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["E",14],Style["σ(E)",14]},PlotLegends->Placed[{Style["",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]

Table[NIntegrate[ψi[x]*ψg[x,t],{x,-∞,∞}],{t,1,2}]

wavecorr={0.8571209590498464` -0.5002658633645123` I,0.4776514950785147` -0.8443484831610817` I,-0.01076557475377125`-0.9339626674721037` I,-0.45235761313091666`-0.7616931731745453` I,-0.7212171220838345`-0.4063642542290182` I,-0.7622384760803487`+0.004148082691040864` I,-0.6002711831761578`+0.3438026701238057` I,-0.3184067053792045`+0.530685128785455` I,-0.019216187467171887`+0.5456335801642568` I,0.21279201797736072` +0.4247743352499456` I}

ListPlot[(Tooltip[{Re[#1],Im[#1]}]&)/@wavecorr,AspectRatio->1]

Inner[5,5,x]

ListPlot[(Tooltip[{Re[#1],Im[#1]}]&)/@%171,AspectRatio->1]

pplus=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,d,100}]},{t,1,50,2}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,d,100000}]}],
{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,d,100000}]}]
p1=ListLinePlot[%,PlotStyle->{Thick},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(>)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
pdwell=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,-d,d}]},{t,1,50,5}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,-d,d}]}],{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,-d,d}]}]
p2=ListPlot[%,PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(<>)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
pminus=Append[Append[Table[{t,NIntegrate[Abs[ψg[x,t]]^2,{x,-100,-d}]},{t,1,50,2}],{100,NIntegrate[Abs[ψg[x,100]]^2,{x,-100000,-d}]}],{1000,NIntegrate[Abs[ψg[x,1000]]^2,{x,-100000,-d}]}]
p3=ListLinePlot[%,PlotStyle->{Dashed},PlotTheme->"Scientific",FrameLabel->{Style["t (ps)",14],Style["",14]},PlotLegends->Placed[{Style["P(<)",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Show[p1,p2,p3]

ClearAll["Global`*"];
β=1000;
S[z_]=Sqrt[((β/π))]Exp[-β* z^2];
Integrate[S[z],{z,-Infinity,Infinity}]
Plot[S[z],{z,-6,6},PlotRange->Full];
NIntegrate[S[x]*Cos[x],{x,-10,10}]
ClearAll[k1,k2,k3,σ];
σ=.02;(*Gaussian params*)
z0=8;
ψ[z_]=((1/(2 π σ^2)))^.25 Exp[((-(z+z0)^2)/(4 σ^2))];
ϵ=0.03;(*distance between consecutive deltas*)
b=-z0;(*central delta position given by mean of Gaussian*)
a=-ϵ+b;
c=ϵ+b;
Plot[ψ[z],{z,-8,8},PlotRange->All];
Integrate[ψ[z]^2,{z,-Infinity,Infinity}]
k3=NIntegrate[ψ[z],{z,-11,a}];(*(8 π σ^2)^.5*)
k1=NIntegrate[ψ[z],{z,a,c}]
k2=NIntegrate[ψ[z],{z,c,10}]
P[z_,0]=k3*S[z-a]+k1*S[z-b]+k2*S[z-c];
Plot[{P[z,0],ψ[z]},{z,-z0-.2,-z0+.2},PlotRange->All]
Plot[{{k3*S[z-a]},{k1*S[z-b]},{k2*S[z-c]},{Exp[-(z-0)^2/(4*σ^2)]/(4*Pi*σ^2)}},{z,-1.4,-0.6},PlotRange->All];
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
Clear[σ];
σ=0;
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
d=0.5;
k0=3;
ν0=0;

k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
f1[x_,t_]=(m/(ℏ*t))(Abs[x]*UnitStep[Abs[x]-d]+d*UnitStep[d-Abs[x]]+xi-2d);
ψg[x_,xi_,t_]=UnitStep[Abs[x]-d]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(x-UnitStep[x-d]2d+xi)^2)/(2(ℏ*t)))]+UnitStep[d-Abs[x]]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(-d+xi)^2)/(2(ℏ*t)))]+((1/(2π)))*((I*m*2*d)/(ℏ^2))*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f1[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f1[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+k[1])])+I*π*Exp[(-I f1[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+k[-1])])))];
t=1;
V0=10;
p1=Plot[{Abs[k3*ψg[x,a,t]+k2*ψg[x,b,t]+k1*ψg[x,c,t]]^2},{x,-8,8},PlotRange->Full,PlotStyle->{RGBColor[0.5,0,0.5],Thick}]

V0=10;
Print["The pole value is"]
k[1]
Plot3D[Abs[ψg[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψg[x,t]]^2,{x,-30,30},{t,0,50}]
Plot[(1/(1^2+(((m*ℏ (En-V0)2d)/(ℏ^2 *Sqrt[2 m En])))^2)),{En,0,30},PlotLegends->"T(E)"]
t=1;
p1=Plot[{Abs[ψg[x,t]]^2},{x,-100,100},PlotRange->Full,PlotStyle->{RGBColor[0.5,0,0.5],Thick}]
Clear[t]
t=25;
p2=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotRange->Full,PlotStyle->{RGBColor[0,0,1],Dotted}]
Clear[t]
t=100;
p3=Plot[{Abs[ψg[x,t]]^2},{x,-100,100},PlotRange->Full,PlotStyle->{RGBColor[1,0,0],Dashed}]
plots=Show[p1,p2,p3]
Style[plots,ImageSize->1]
Clear[V0,t,k];
Remove[k[i]]
k[i_]:=((1+i*1*Sqrt[1+8d^2m*V0/ℏ^2])/(2 d));
V0=10;
t=25;
Print["The pole value is"]
k[1]
p4=Plot[{Abs[ψg[x,t]]^2},{x,-10,20},PlotRange->Full,PlotStyle->{RGBColor[0,0,1],Dotted}]
Show[p4,p2]
Clear[t]
NIntegrate[Abs[ψg[x,1]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,10]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,50]]^2,{x,-100,100},MaxRecursion->20]

Let us introduce the wavefunction of the delta-function intial wave packet as ψd

ClearAll["Global`*"]
Remove["Global`*"]
Unprotect["Global`*"]
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=0;
ℏ=1.055;
m=9.109;
k0=3;
xi=8;
d=0.5;
k0=3;
ν0=0;

k[i_]:=((1+i*1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
f[x_,t_]=(1/(2(σ^2+((I*ℏ*t)/(2m)))))(I*xi+I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
f1[x_,t_]=(m/(ℏ*t))(Abs[x]*UnitStep[Abs[x]-d]+d*UnitStep[d-Abs[x]]+xi-2d);
ψg[x_,xi_,t_]=UnitStep[Abs[x]-d]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(x-UnitStep[x-d]2d+xi)^2)/(2(ℏ*t)))]+UnitStep[d-Abs[x]]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(-d+xi)^2)/(2(ℏ*t)))]+((1/(2π)))*((I*m*2*d)/(ℏ^2))*op[t][(Exp[(σ^2+((I*ℏ*t)/(2m)))(f1[x,t])^2]*((1/(k[1]-k[-1])))(-I*π*Exp[(-I f1[x,t]+k[1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+k[1])])+I*π*Exp[(-I f1[x,t]+k[-1])^2((σ^2)+((I*ℏ*t)/(2m)))](Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+k[-1])])))];
t=1;
V0=10;
p1=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotRange->Full,PlotStyle->{RGBColor[0.5,0,0.5],Thick}]

V0=10;
Print["The pole value is"]
k[1]
Plot3D[Abs[ψg[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψg[x,t]]^2,{x,-30,30},{t,0,50}]
Plot[(1/(1^2+(((m*ℏ (En-V0)2d)/(ℏ^2 *Sqrt[2 m En])))^2)),{En,0,30},PlotLegends->"T(E)"]
t=1;
p1=Plot[{Abs[ψg[x,t]]^2},{x,-100,100},PlotRange->Full,PlotStyle->{RGBColor[0.5,0,0.5],Thick}]
Clear[t]
t=25;
p2=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotRange->Full,PlotStyle->{RGBColor[0,0,1],Dotted}]
Clear[t]
t=100;
p3=Plot[{Abs[ψg[x,t]]^2},{x,-100,100},PlotRange->Full,PlotStyle->{RGBColor[1,0,0],Dashed}]
plots=Show[p1,p2,p3]
Style[plots,ImageSize->1]
Clear[V0,t,k];
Remove[k[i]]
k[i_]:=((1+i*1*Sqrt[1+8d^2m*V0/ℏ^2])/(2 d));
V0=10;
t=25;
Print["The pole value is"]
k[1]
p4=Plot[{Abs[ψg[x,t]]^2},{x,-10,20},PlotRange->Full,PlotStyle->{RGBColor[0,0,1],Dotted}]
Show[p4,p2]
Clear[t]
NIntegrate[Abs[ψg[x,1]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,10]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,50]]^2,{x,-100,100},MaxRecursion->20]

Some useful integrals

(1/(Sqrt[2 π]))({{\[Piecewise], ({{((2 (c^(3/2)))/((c-I k) (c+I k))), Im[k]+Re[c]>0&&Im[k]-Re[c]<0}, {Integrate[(Sqrt[c]) (E^(-I k x-c Abs[x])),{x,-∞,∞},Assumptions->Im[k]-Re[c]≥0||Im[k]+Re[c]≤0], True}})}})

ClearAll["Global`*"];
a=I ℏ t/2m;
Integrate[Exp[-a u^2]/(u+ c+I b),{u,-∞,∞},Assumptions->b>0&&m ∈Reals &&ℏ ∈Reals&&t ∈Reals]
Integrate[Exp[-a u^2]/(u+ c-I b),{u,-∞,∞},Assumptions->b>0&&m ∈Reals &&ℏ ∈Reals&&t ∈Reals]
Integrate[Exp[-a u^2]/(u+ c+I b),{u,-∞,∞},Assumptions->b<0&&m ∈Reals &&ℏ ∈Reals&&t ∈Reals]
Integrate[Exp[-a u^2]/(u+ c-I b),{u,-∞,∞},Assumptions->b<0&&m ∈Reals &&ℏ ∈Reals&&t ∈Reals]
Integrate[Exp[-a u^2],{u,-∞,∞},Assumptions->b>0&&c>0&&m ∈Reals &&ℏ ∈Reals&&t ∈Reals]

What follows is the time evolution of an initial wave-packet is exponential;

Remove["Global`*"];
Unprotect["Global`*"];
ClearAll["Global`*"];
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
σ=0;
f1[x_,t_]=(1/(2(((I*ℏ*t)/(2m)))))(I(Abs[x]*HeavisideTheta[Abs[x]-d]+d*UnitStep[d-Abs[x]]-2d));
(Subscript[c, 1])=-I((1+1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
(Subscript[c, 2])=-I((1-1*Sqrt[1-8d^2m*V0/ℏ^2])/(2 d));
(Subscript[c, 3])=I c0;
(Subscript[c, 4])=-I c0;
ψg[x_,t_]=((UnitStep[RealAbs[x]-d])/π)c0^{3/2}((Exp[-(( m)/(I 2ℏ t))(x-UnitStep[x-d]2d)^2])/(2I c0))(-I π Exp[((I ℏ t)/(2m))( -c0- I f1[x,t])^2]Erfc[Sqrt[((I ℏ t)/(2m))](- c0-  I f1[x,t])]+I π Exp[((I ℏ t)/(2m))(c0-I f1[x,t])^2]Erfc[Sqrt[((I ℏ t)/(2m))]( c0-I f1[x,t])]);
c0=((m V0 2 d)/(ℏ^2));
m=9.109;
ℏ=1.055;
d=0.5;
V0=10;
t=5;
ψg[5,t]
(*Plot3D[Abs[ψg[x,t]]^2,{x,-20,20},{t,0,50}]
ContourPlot[Abs[ψg[x,t]]^2,{x,-30,30},{t,0,50}]
Plot[(1/(1^2+(((m*ℏ (En-V0)2d)/(ℏ^2 *Sqrt[2 m En])))^2)),{En,0,30},PlotLegends->"T(E)"]*)

p1=Plot[{Abs[ψg[x,t]]^2},{x,-10,8},PlotStyle->{RGBColor[0,0,1],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=5",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+1}}]]
(*Clear[t]
t=10
p2a=Plot[{Abs[ψg[x,t]]^2},{x,-8,8},PlotStyle->{RGBColor[0,0,0],Dotted},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=10",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5}}]]
Clear[t]
t=25;
p3=Plot[{Abs[ψg[x,t]]^2},{x,-20,20},PlotStyle->{RGBColor[0,0,1],Dashed},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=25",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5+.5}}]]
Clear[t]
t=50;
p4=Plot[{Abs[ψg[x,t]]^2},{x,-20,20},PlotStyle->{RGBColor[1,1,0],Thick},PlotTheme->"Scientific",PlotRange->Full,PlotLabel->{Style["Probability distribution versus position",14]},FrameLabel->{Style["x (nm)",14],Style["|ψ(x,t)((|)^2)",14]},PlotLegends->Placed[{Style["t=50",13]},{Scaled[{0,0.5}],{-1.4+0.2,-.7-.5-.5}}]]
plots=Show[p3,p2a,p4,p1]
Style[plots,ImageSize->1]
NIntegrate[Abs[ψg[x,1]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,10]]^2,{x,-100,100},MaxRecursion->20]
NIntegrate[Abs[ψg[x,50]]^2,{x,-100,100},MaxRecursion->20]*)

+((UnitStep[d-RealAbs[x]])/π)c0^{3/2}((Exp[-(( m)/(I 2ℏ t))(d)^2])/(2I c0))(-I π Exp[((I ℏ t)/(2m))( -c0-I f1[x,t])^2]Erfc[Sqrt[((I ℏ t)/(2m))](- c0-I f1[x,t])]+I π Exp[((I ℏ t)/(2m))(c0-I f1[x,t])^2]Erfc[Sqrt[((I ℏ t)/(2m))]( c0-I f1[x,t])])+0((c0^(3/2))/π)(((I*m*2*d)/(ℏ^2)))op[t][Exp[((I*ℏ*t)/(2m))(f1[x,t])^2]((1/(((Subscript[c, 1])-(Subscript[c, 2]))((Subscript[c, 1])-(Subscript[c, 3]))((Subscript[c, 1])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 1]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 1]))]+(1/(((Subscript[c, 2])-(Subscript[c, 1]))((Subscript[c, 2])-(Subscript[c, 3]))((Subscript[c, 2])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 2]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 2]))]+(1/(((Subscript[c, 3])-(Subscript[c, 1]))((Subscript[c, 3])-(Subscript[c, 2]))((Subscript[c, 3])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 3]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 3]))]+(1/(((Subscript[c, 4])-(Subscript[c, 1]))((Subscript[c, 4])-(Subscript[c, 2]))((Subscript[c, 4])-(Subscript[c, 3]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 4]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 4]))])];

Erfc[Sqrt[((I ℏ t)/(2m))](c0-I f1[x,t])]/.x->10

s

+((UnitStep[-RealAbs[x]+d])/π)c0^{3/2}Exp[-((I ℏ)/(2m t))(((m d)/(ℏ t)))^2](-I π Exp[((I ℏ t)/(2m))(c0-((I m d)/(ℏ t)))^2])+0((c0^(3/2))/π)(((I*m*2*d)/(ℏ^2)))op[t][Exp[((I*ℏ*t)/(2m))(f1[x,t])^2](((0*1)/(((Subscript[c, 1])-(Subscript[c, 2]))((Subscript[c, 1])-(Subscript[c, 3]))((Subscript[c, 1])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]-I (Subscript[c, 1]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]-I (Subscript[c, 1]))]+((0*1)/(((Subscript[c, 2])-(Subscript[c, 1]))((Subscript[c, 2])-(Subscript[c, 3]))((Subscript[c, 2])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]-I (Subscript[c, 2]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]-I (Subscript[c, 2]))]+(1/(((Subscript[c, 3])-(Subscript[c, 1]))((Subscript[c, 3])-(Subscript[c, 2]))((Subscript[c, 3])-(Subscript[c, 4]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 3]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 3]))]+(1/(((Subscript[c, 4])-(Subscript[c, 1]))((Subscript[c, 4])-(Subscript[c, 2]))((Subscript[c, 4])-(Subscript[c, 3]))))(-I π Exp[(-I f1[x,t]+I (Subscript[c, 4]))^2((σ^2)+((I*ℏ*t)/(2m)))])Erfc[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I f1[x,t]+I (Subscript[c, 4]))])];
c0=((m V0 2 d)/(ℏ^2));
m=9.109;
ℏ=1.055;
d=0.5;
V0=10;
t=5;
ψg[5,t]

ψg[x,t]

(-I π Exp[((I ℏ t)/(2m))(c0-((I m d)/(ℏ t)))^2])+((c0^(3/2))/π)(((I*m*2*d)/(ℏ^2)))op[t][Exp[((I*ℏ*t)/(2m))(f1[x,t])^2]Sum[Product[((KroneckerDelta[n,j]((Subscript[c, n])-(Subscript[c, n])))/(((Subscript[c, n])-(Subscript[c, j])))),{j,1,4}](-I*π*Exp[(I*f[x,t]-I*(Subscript[c, i]))^2Sqrt[((I*ℏ*t)/(2m))]](-Sign[Im[((Subscript[c, i])-I*f[x,t])]]+Erf[Sqrt[((I*ℏ*t)/(2m))](-I(f[x,t]-(Subscript[c, i])))])),{n,1,4}]]

ClearAll["Global`*"];
Factor[-I d x^2+x+I ((m v0)/(ℏ^2))2d]

Solve[-I d x^2+x+I ((m v0)/(ℏ^2))2d==0,x]

kosuru calculations;

FourierTransform[Sqrt[c]*Exp[-c*Abs[x]],x,k]
1/Sqrt[2π]*Integrate[Sqrt[c]*Exp[-c*Abs[x]-I*k*x],{x,-∞,∞}]

Clear[d];
d=2;
Plot[Abs[x]*UnitStep[Abs[x]-d]+d *UnitStep[-Abs[x]+d],{x,-3,3}]

Deriving the kernels

ClearAll["Global`*"];
y1[x_]=Exp[I*k*x]+b*Exp[-I*k*x];(*Asymptotic wave before lens1*)
y2[x_]=ϕ;
y5[x_]=c*Exp[I*k*x];(*After lens2*)
Solve[{y1[-d]==y2[-d],y2[d]==y5[d],y5'[d]-y1'[-d]==-2(m/(ℏ^2))(En-V)2d*y2[0]},{b,c,ϕ}]
{{b,c,ϕ}}={b,c,ϕ}/.%

f1[x_,t_]=(m/(ℏ*t))(Abs[x]*UnitStep[Abs[x]-d]+d*UnitStep[d-Abs[x]]+xi-2d);
ψd[x_,t_]=UnitStep[Abs[x]-d]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(x-UnitStep[x-d]2d+xi)^2)/(2(ℏ*t)))]+UnitStep[d-Abs[x]]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(d+xi)^2)/(2(ℏ*t)))]+Sign[x+d]((1/(2π)))*((I*m*2*d)/(ℏ^2))*op[t][Exp[(((I*ℏ*t)/(2m)))(f1[x,t])^2]((1/((Subscript[k, 1])-(Subscript[k, -1]))))*Sum[(2i-1)*I*π*Exp[(I*f1[x,t]-I*(Subscript[k, 2i-1]))^2Sqrt[((I*ℏ*t)/(2m))]](-Sign[Im[(Subscript[k, 2i-1])+f1[x,t]]]+Erf[Sqrt[((I*ℏ*t)/(2m))](-I(f1[x,t]+(Subscript[k, 2i-1])))]),{i,0,1}]];
ψg[x_,t_]=UnitStep[Abs[x]-d](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi)-((ℏ k0^2t)/(2m)))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+UnitStep[d-Abs[x]](((σ^2)/(2π(σ^2+((I*ℏ*t)/(2m)))^2)))^0.25*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi)-((ℏ k0^2t)/(2m)))-(((d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+Sign[x+d](((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2+I*ν0/ℏ]*op[t][Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]((1/((Subscript[k, 1])-(Subscript[k, -1]))))*Sum[-(2i-1)*I*π*Exp[(I*f[x,t]-I*(Subscript[k, 2i-1]))^2((σ^2)+((I*ℏ*t)/(2m)))](Sign[κ[x,t]+Im[(Subscript[k, 2i-1])]]-Erf[Sqrt[σ^2+((I*ℏ*t)/(2m))](κ[x,t]+Im[(Subscript[k, 2i-1])])]),{i,0,1}]]

ComplexExpand[1/(a+I b)]

ψg[x,5]

What follows is the time evolution of an initial wave-packet is exponential;

ClearAll["Global`*"];
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
ψe[x_,t_]=((c0^(3/2))/π)Sign[x+d](((I*m*2*d)/(ℏ^2)))Exp[((I*ℏ*t)/(2m))(f1[x,t])^2]op[t][Sum[Product[((KroneckerDelta[n,j]((Subscript[c, n])-(Subscript[c, n])))/(((Subscript[c, n])-(Subscript[c, j])))),{j,1,4}](-I*π*Exp[(I*f[x,t]-I*(Subscript[c, i]))^2Sqrt[((I*ℏ*t)/(2m))]](-Sign[Im[((Subscript[c, i])-I*f[x,t])]]+Erf[Sqrt[((I*ℏ*t)/(2m))](-I(f[x,t]-(Subscript[c, i])))])),{n,1,4}]]

Simplify[Integrate[Exp[-a*x^2]/(x+I*c+f),{x,-∞,∞},Assumptions->Re[a]>0&&Element[c,Reals]&&c>0&&Element[f,Reals]]]

PLot of the |Green’s function((|)^2) (initial δ(x+1)) versus x

ClearAll["Global`*"];
Remove["Global`*"];
Unprotect["Global`*"];
m=1;
ℏ=1;
V0=1;
d=.7;
xi=a;
a=3;
(Subscript[k, i_])=((-I+i*I*Sqrt[1+8d^2m*V0/ℏ^2])/(2d));

op[t_]=(I*ℏ*(∂_t)#-V0#)&;
f1[x_,t_]=(m/(ℏ*t))(Abs[x]*UnitStep[Abs[x]-d]+d*UnitStep[d-Abs[x]]-a-2d);
f2[x_,t_]=Simplify[Exp[(((I*ℏ*t)/(2m)))(f1[x,t])^2]((1/((Subscript[k, 1])-(Subscript[k, -1]))))*Sum[(2i-1)*I*π*Exp[(I*f1[x,t]-I*(Subscript[k, 2i-1]))^2Sqrt[((I*ℏ*t)/(2m))]](-Sign[Im[(Subscript[k, 2i-1])+f1[x,t]]]+Erf[Sqrt[((I*ℏ*t)/(2m))](-I(f1[x,t]+(Subscript[k, 2i-1])))]),{i,0,1}],Assumptions->Element[x,Reals]&&Element[t,Reals]&&t>0];
G[x_,a_,t_]=UnitStep[Abs[x]-d]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(x-UnitStep[x-d]2d-xi)^2)/(2(ℏ*t)))]+UnitStep[d-Abs[x]]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(d+xi)^2)/(2(ℏ*t)))]+Sign[x+d]((1/(2π)))*((I*m*2*d)/(ℏ^2))*op[t][f2[x,t]];
(*Plot[{Re[G[x,a,t]],Im[G[x,a,t]]},{x,-40,40},PlotRange->{-10,10}];
ψsq[x_,t_]=Evaluate[Abs[G[x,a,t]]^2,Element[x,Reals]];Assumptions->Element[m,Reals]&&m>0&&Element[ℏ,Reals]&&ℏ>0&&Element[d,Reals]&&d>0&&Element[xi,Reals]&&Element[a,Reals]&&Element[V0,Reals]];*)
xp=-1.5;
yp=-.5;
ls=12;
sh=0;
t=2;
Plot[{Abs[G[x,a,t]]^2},{x,-100,100},PlotRange->All,PlotStyle->{Thick},PlotLegends->Placed[{"|G[x,3,2]((|)^2)"},{Scaled[{0,0.5}],{yp-1,xp-2sh+.5}}],PlotTheme->{"Scientific","GrayColor"},FrameLabel->{Style["time (ps)",ls],Style["Probability amplitude |ψ|^2",ls]},AxesStyle->Thick]
NIntegrate[Abs[G[x,a,t]]^2,{x,-1000,1000}]

PLot of the |ψ(x,t)((|)^2) for the initial Gaussian wavepacket

ClearAll["Global`*"];
Remove["Global`*"]
Unprotect["Global`*"];
m=1;
ℏ=1;
σ=0.5;
V0=1;
d=.5;
k0=1;
ν0=1;
xi=6;
(Subscript[k, i_])=((-I+i*I*Sqrt[1+8d^2m*V0/ℏ^2])/(2d));
op[t_]=(I*ℏ*(∂_t)#-V0#)&;
f[x_,t_]=(((σ^2-((I*ℏ*t)/(2m))))/(2(σ^4+(((ℏ*t)/(2m)))^2)))(2k0*(σ^2)+I*xi);
Table[Simplify[Sign[Abs[x]*UnitStep[Abs[x]-d]+Im[(Subscript[k, 1])+f[x,t]]],Assumptions->Element[t,Reals]&&t>0&&t<2.1623050908562145`],{x,-50,50}]
Table[Simplify[Sign[Abs[x]*UnitStep[Abs[x]-d]+Im[(Subscript[k, -1])+f[x,t]]],Assumptions->Element[t,Reals]&&t>0&&t<2.1623050908562145`],{x,-50,50}]
(*Table[Simplify[Sign[Im[(Subscript[k, -1])+f[x,t]]],Assumptions->Element[t,Reals]&&t>0],{x,-50,50}]*)
Clear[x,t];
f2[x_,t_]=Simplify[Sum[(2i-1)*I*π*Exp[(I*f[x,t]-I*(Subscript[k, 2i-1]))^2Sqrt[(σ^2)+((I*ℏ*t)/(2m))]](-1*KroneckerDelta[i,1]+Sign[Abs[x]-2.1623050908562145`]*KroneckerDelta[i,0]+Erf[Sqrt[σ^2+((I*ℏ*t)/(2m))](-I(f[x,t]+(Subscript[k, 2i-1])))]),{i,0,1}],Assumptions->Element[x,Reals]&&Element[t,Reals]&&t>0];
ψg[x_,t_]=UnitStep[Abs[x]-d]((σ/((2π)^0.5(σ^2+((I*ℏ*t)/(2m))))))^0.5*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi))-((ℏ k0^2t)/(2m))-(((x-UnitStep[x-d]2d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+UnitStep[d-Abs[x]]((σ/((2π)^0.5(σ^2+((I*ℏ*t)/(2m))))))^0.5*Exp[I*((σ^2)/(σ^2+((I*ℏ*t)/(2m))))(k0(x+xi))-((ℏ k0^2t)/(2m))-(((d+xi)^2)/(4(σ^2+((I*ℏ*t)/(2m)))))+((I*ν0)/ℏ)]+Sign[x+d](((σ^2)/(2(π^3))))^.25*((I*m*2*d)/(ℏ^2))*Exp[-σ^2k0^2+I*ν0/ℏ]*op[t][((1/((Subscript[k, 1])-(Subscript[k, -1]))))*Exp[(σ^2+((I*ℏ*t)/(2m)))(f[x,t])^2]* f2[x,t]]
Print["Check with initial condition ψ(t)-ψ(0)"]
Simplify[ψg[x,0]-((1/(Sqrt[2*π*σ^2])))^.5*Exp[-(((x+xi)^2)/(4σ^2))+I*k0(x+xi)+I*((ν0)/ℏ)]]
Clear[x];
t=0;
Plot[{Abs[ψg[x,t]]^2},{x,-40,40},PlotRange->Full,Exclusions->None]
p1=Plot[{Abs[((1/(Sqrt[2*π*σ^2])))^.5*Exp[-(((x+xi)^2)/(4σ^2))+I*k0(x+xi)+I*((ν0)/ℏ)]]^2},{x,-40,40},PlotRange->Full,Exclusions->None]
p2=Plot[{Abs[((1/(Sqrt[2*π*σ^2])))^.5*UnitStep[Abs[x]-d]*Exp[-(((x-UnitStep[x-d]*2d+xi)^2)/(4σ^2))+I*k0(x+xi)+I*((ν0)/ℏ)]]^2},{x,-40,40},PlotRange->Full,Exclusions->None]
Show[p1,p2]

Table[{Sign[Abs[x]-3]},{x,-40,40}]

Numerical evaluation of the wavepacket using ND solve of mthematica

Unprotect["Global'*"];
ClearAll["Global'*"];
Remove["Global`*"];
x0=1;
σ=1/Sqrt[40];
m=1;
ℏ=1;
V0=1;
d=-.7;
dThick=40;
dTime=100;
eqSimplee=I*ℏ*D[Psi[x,t],t]==-((ℏ^2)/(2m))*D[Psi[x,t],{x,2}]+V0*UnitStep[x+d]*UnitStep[d-x]*Psi[x,t];
β=1/(4σ^2);
eqInie=Psi[x,0]==Sqrt[(β/π)]*Exp[-β*(x+x0)^2];
eqBlefte=Re[Psi[-dThick,t]]==0+0*I;
eqBrighte=Re[Psi[dThick,t]]==0+0*I;
h=-.4;
y=-1.5;
ys=0;
(*mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};*)
res=NDSolve[{eqSimplee,eqInie,eqBlefte,eqBrighte},Psi[x,t],{x,-dThick,dThick},{t,0,dTime}]
Plot3D[Re[Psi[x,t]/.res],{x,-dThick,dThick},{t,0,dTime},PlotRange->Full]
t=0;
aa1=Plot[Evaluate[Abs[Psi[x,t]/.res]^2],{x,-dThick,dThick},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, g])(t)",14]},PlotStyle->{Dashed}]

Plot3D[Re[Psi/.%],{x,-dThick,dThick},{t,0,dTime}];

ClearAll["Global`*"];
(Subscript[k, j_])=((-I+j*I*Sqrt[1+8d^2m*V0/ℏ^2])/(2d));
Needs["NumericalCalculus`"]

op[t_]=(I*ℏ*ND[#,t,2]-V0#)&;
f1[x_,t_]=(m/(ℏ*t))(Abs[x]*UnitStep[Abs[x]-d]+d*UnitStep[d-Abs[x]]+xi-2d);

Ganalytica[x_,s_]=(((-m)/(I*ℏ^2)))Sqrt[(ℏ/(2m*I*s))]Exp[I*Sqrt[((2m*I*s)/ℏ)]Abs[x-x0]]+(1/(1+(m/(I*ℏ^2))(I*ℏ*s+V0)2d*Sqrt[(ℏ/(2m*I*s))]Exp[-I*Sqrt[((2m*I*s)/ℏ)]Abs[d]]))(I*ℏ*s+V0)2d(((((-m)/(I*ℏ^2))*Sqrt[(ℏ/(2m*I*s))]))^2)Exp[I*Sqrt[((2m*I*s)/ℏ)](Abs[x]+Abs[x0]-Abs[d])];
G1[x_,s_]=LaplaceTransform[UnitStep[Abs[x]-d]*((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(x-UnitStep[x-d]2d+xi)^2)/(2(ℏ*t)))],t,s]
G2[x_,s_]=LaplaceTransform[UnitStep[d-Abs[x]]((m/(2π(I*ℏ*t))))^0.5*Exp[-((I*m(d-xi)^2)/(2(ℏ*t)))],t,s]
m=1;
ℏ=1;
V0=1;
d=1;
x0=-xi;
xi=1;
x=1;
a=1;
p2=Plot[Re[Ganalytica[x,s]],{s,0,30},PlotRange->All]
g3[x_,t_]=Sign[x+d]((1/(2π)))*((I*m*2*d)/(ℏ^2))*op[t][Exp[(((I*ℏ*t)/(2m)))(f1[x,t])^2]((1/((Subscript[k, 1])-(Subscript[k, -1]))))*Sum[(2i-1)*I*π*Exp[(I*f1[x,t]-I*(Subscript[k, 2i-1]))^2((I*ℏ*t)/(2m))](-Sign[Im[(Subscript[k, 2i-1])+f1[x,t]]]+Erf[Sqrt[((I*ℏ*t)/(2m))](-I(f1[x,t]+(Subscript[k, 2i-1])))]),{i,0,1}]];
Table[{s,Re[NIntegrate[g3[x,t]*Exp[-s*t],{t,1,100}]+G1[x,s]+G2[x,s]]},{s,0,30}]
p1=ListPlot[%,PlotJoined->True]
Show[p1,p2]

Checking the coorectness of the Laplace domain green's function given in deepak's paper;

Ganalytica[x_,s_]=(((-m)/(I*ℏ^2)))Sqrt[(ℏ/(2m*I*s))]Exp[I*Sqrt[((2m*I*s)/ℏ)]Abs[x-x0]]+(1/(1+(m/(I*ℏ^2))(I*ℏ*s-V0)2d*Sqrt[(ℏ/(2m*I*s))]Exp[-I*Sqrt[((2m*I*s)/ℏ)]Abs[d]]))(I*ℏ*s-V0)2d(((((-m)/(I*ℏ^2))*Sqrt[(ℏ/(2m*I*s))]))^2)Exp[I*Sqrt[((2m*I*s)/ℏ)](Abs[x]+Abs[x0]-Abs[d])];

Apart[1/((k-c1)(k-c2)(k-c3)(k-c4))]

ClearAll["Global`*"];
Integrate[Exp[-a*x^2]/(x+I*b+c),{x,-∞,∞},Assumptions->Element[b,Reals]&&Element[c,Reals]&&b>0]
Integrate[Exp[-a*x^2]/(x+I*b+c),{x,-∞,∞},Assumptions->Element[b,Reals]&&b<0]

ClearAll["Global`*"]
V0=1;
d=1;
a=1;
Plot[{HeavisideTheta[x+d]*V0*HeavisideTheta[d-x],Sqrt[a/π]*Exp[-a*(x+7)^2]},{x,-d-10,d+2},PlotRange->Full,Exclusions->None,Ticks->None,Axes->{True,False}]

ClearAll["Global`*"];
y1[x_]=A1*Exp[Sqrt[s/d]*x];
y2[x_]=A2*Exp[-Sqrt[s/d]*x]+A3*Exp[+Sqrt[s/d]*x];
y3[x_]=A4*Exp[-Sqrt[s/d]*x];
Solve[{y1[x1]==y2[x1],y2'[x1]-y1'[x1]==-k1/d,y2[xc]==y3[xc],y3'[xc]-y2'[xc]==-k1/d*y2[xc]},{A1,A2,A3,A4}]

ClearAll["Global`*"];
s=FourierTransform[Sqrt[c0]*Exp[-c0*Abs[x]],x,k]
c0=2;
k=2;
Integrate[1/Sqrt[2π]Sqrt[c0]*Exp[-c0*Abs[x]-I*k*x],{x,-∞,∞}]
s

ClearAll["Global`*"];
Integrate[(((c0^(3/2)) (Sqrt[(2/π)]))/((c0^2)+(k^2)))UnitStep[Abs[x]-d]Exp[-I*((ℏ*k^2)/(2m))t+I*k*(x-a)],{k,-∞,∞}]

ClearAll["Global`*"];
V0=0;
ℏ=1;
m=1;
d=1;
Plot[(1/(1+(((m(En-V0)2 d)/(ℏ^2 Sqrt[2m En]/ℏ)))^2)),{En,0,10},PlotRange->All]