(* here are the dimensionless parameters of our wave packet *)
(* inverse width of the wave packet in momentum space *) a=0.5;
(* parameter k0 *) k0=1;
(* cutoff for momentum integrals *) cutoff=4;
(* initial position of the particle is centered around x0 *) x0=-20;
(* parameter Λ *) Λ=1;
(* K can be expressed in terms of k and Λ *) κ[k_]:=Sqrt[k^2-Λ^2]

(* now specify the relevant energy eigenstates *)
(* transmission amplitude *)
tra[k_]:=2*κ[k]*k*Exp[-2*I*k]/(-I*(k^2+κ[k]^2)*Sin[2*κ[k]]+2*k*κ[k]*Cos[2*κ[k]])
(* reflection amplitude *)
ref[k_]:=-I(k^2-κ[k]^2)*Exp[-2*I*k]*Sin[2*κ[k]]/(-I*(k^2+κ[k]^2)*Sin[2*κ[k]]+2*k*κ[k]*Cos[2*κ[k]])

(* amplitude A*) 
aa[k_]:=Exp[-I*κ[k]]/2/κ[k]*(Exp[-I*k]*(κ[k]-k)+ref[k]*Exp[I*k]*(κ[k]+k));

(* amplitude B *)
bb[k_]:=Exp[I*κ[k]]/2/κ[k]*(Exp[-I*k]*(κ[k]+k)+ref[k]*Exp[I*k]*(κ[k]-k));

(* un-normalized wave function of scattering state in the three different spatial regions *)

(* x<-a*) 
psi1[x_,t_]:=NIntegrate[Exp[-a*(k-k0)^2-I*k*x0]*(Exp[I*k*x]+ref[k]*Exp[-I*k*x])*Exp[-I*t*k^2/2],{k,Λ,cutoff}];

(* -a<x<a *)
psi2[x_,t_]:=NIntegrate[Exp[-a*(k-k0)^2-I*k*x0]*(aa[k]*Exp[-I*κ[k]*x]+bb[k]*Exp[I*κ[k]*x])*Exp[-I*t*k^2/2],{k,Λ,cutoff}];

(* -x>a *)
psi3[x_,t_]:=NIntegrate[Exp[-a*(k-k0)^2-I*k*x0]*tra[k]*Exp[I*k*x]*Exp[-I*t*k^2/2],{k,Λ,cutoff}];

(* dimensionless unnormalized probability density *)
prob0[x_,t_]:=If[x<=-1,Abs[psi1[x,t]]^2,If[x<=1,Abs[psi2[x,t]]^2,Abs[psi3[x,t]]^2]]

(* calculate it at time tau=0 *)
data=Table[{0.125*j,prob0[0.125*j,0]},{j,-240,0}];
(* interpolation function *)
px0=Interpolation[data];
(* calculate the square of normalization factor *)
nc=NIntegrate[px0[x],{x,-30,-10}];

(* dimensionless normalized probability density *)prob[x_,t_]:=prob0[x,t]/nc;

(* Plotting the solution*)
(* first look at the transmission and reflection amplitudes as functions of wave number*)
Plot[{Abs[ref[k]]^2,Abs[tra[k]]^2},{k,Λ,cutoff}]

(* now plot the probability density of the initial wave packet *)

data=Table[{0.125*j,prob[0.125*j,0]},{j,-240,0}];
ListPlot[data,PlotRange->All]

(* check that this is indeed normalized *)
px0=Interpolation[data];
NIntegrate[px0[x],{x,-30,-10}]

data=Table[{0.125*j,prob[0.125*j,10]},{j,-240,80}];
ListPlot[data,PlotRange->All]

data=Table[{0.15*j,prob[0.15*j,20]},{j,-200,240}];
ListPlot[data,PlotRange->All]

data=Table[{0.175*j,prob[0.175*j,40]},{j,-280,520}];
ListPlot[data,PlotRange->All]

data=Table[{0.25*j,prob[0.25*j,60]},{j,-300,520}];
ListPlot[data,PlotRange->All]