NotebookEvaluate[NotebookDirectory[]<>"/split_step.nb"];

HarmonicTrap1D[];
sizes={4096};
timeStep=0.0001;
k=-0.5;
contactInteractionFactor=4 Pi k;
max=50ho;
min=-50ho;
initialize[];
gridReport

calcAllTable
p1=plotProjections;

AbsoluteTiming[evolve["ite",10,10]]

calcAllTable
p2=plotProjections;
Show[p1,p2,PlotRange->{{-3,3},All}]

Turn off harmonic potential

plot:=AppendTo[lplots,{time,plotProjections,interp[wf]}];
newPotential[(0&)];
reset;
AbsoluteTiming[evolve["ite",10,10,10]]

{ListPlot[lmoments[[All,{1,2}]],Joined->True,PlotRange->All],
ListPlot[lmoments[[All,{1,3}]],Joined->True,PlotRange->All]}

wfsoliton=wf;
psi=interp[wfsoliton]
wfsoliton2=Map[psi[First[#]]&,mesh[],{dim}];
Chop@Norm[wfsoliton-wfsoliton2]

x0=-10;
wfshifted=Map[Piecewise[{{psi[First[#]-x0],Abs[First[#]-x0]<10}},0]&,mesh[],{dim}];

p=3; (* momentum *)
wf=wfshifted*Exp[I p meshX];
ListPlot[projectionX,PlotRange->All,Joined->True]

reset;
AbsoluteTiming[evolve["rte",15,50,10]]

plotM["meanX"]

lmstart=listM["meanX"][[1]]
lmend=listM["meanX"][[-1]]
Δt=lmend[[1]]-lmstart[[1]]
Δx=lmend[[2]]-lmstart[[2]]
v=Δx/Δt

GraphicsGrid[Table[{lplots[[i,2]]},{i,Length[lplots]}],ImageSize->3 72]

(* The slight broadening of the soliton is a numerical artifact due to a finite step size. *)
l=listM["sigmaX"]^2-listM["meanX"]^2;
l[[All,1]]=listM["sigmaX"][[All,1]];
l[[All,2]]=Sqrt[l[[All,2]]];
ListPlot[l,AxesLabel->{"t","(Subscript[σ, x])"}]