Решение уравнения Шрёдингера методом стрельбы

Мы будем решать уравнение Шрёдингера (ℏ = m = 1) в следующем потенциле :

u=3;
U[x_]=-(u/((Cosh[x])^2));

накладывая следующие граничные условия: ψ[x-> -∞]=C (e^(κ x))<-> ψ'[-L]/ψ[-L]=κ, κ=(Sqrt[2|ϵ|]);
L нужно выбрать таким, чтобы потенциал там можно было (численно) выбросить.

L=5;
sol=ParametricNDSolve[{-(1/2)ψ''[x]+U[x] ψ[x]==ϵ ψ[x],
ψ'[-L]==(Sqrt[2Abs[ϵ]]),
ψ[-L]==1},ψ[x],{x,-L,L},ϵ];

Так выглядят решения для различных значений энергии :

Manipulate[Plot[ψ[x][ϵ]/.sol,{x,-(L/2),(L/2)}],{ϵ,-u,0,0.1}]

Plot[Evaluate@Table[ψ[x][ϵ]/.sol,{ϵ,-2.3,-1.7,0.1}],{x,-(L/2),(L/2)},PlotLegends->Table["ϵ="+ϵ,{ϵ,-2.3,-1.7,0.1}]]

Мы хотим, чтобы ψ[x-> +∞]=C1 (e^(-κ x))+C2 (e^(κ x)), и C2=0.  Для этого мы делаем FindFit, который найдёт нам значения C1 и C2

ψAsymptotic[C1_,C2_,ϵ_,x_]=C1 Exp[-(Sqrt[2 Abs[ϵ]])x] + C2 Exp[+(Sqrt[2 Abs[ϵ]])x];
FindAsymptotic=Function[ϵ,
data = Table[{x,ψ[x][ϵ]/.sol},{x,2,L,0.1}];
FindFit[data,ψAsymptotic[C1,C2,ϵ,x],{C1,C2},x]];

Построим рядом функцию и асимптотику, чтобы убедиться, что Fit работает:

With[{ϵ=-0.5},Show[{LogPlot[Abs[ψ[x][ϵ]]/.sol,{x,0,L}],
LogPlot[Evaluate[Abs[ψAsymptotic[C1,C2,ϵ,x]]/.FindAsymptotic[ϵ]],{x,(L/2),L},PlotStyle->Red]}]]

Теперь мы можем строить C2 как функцию ϵ, и искать её нули (соответствующие связанным состояниям):

EigenEquation[ϵ_?NumericQ]:=C2/.FindAsymptotic[ϵ];

Plot[EigenEquation[ϵ],{ϵ,-2.1,0},BaseStyle->{FontSize->20}]

Видно, что имеются два собственных числа:

eig1=FindRoot[EigenEquation[ϵ],{ϵ,-3}]
eig2=FindRoot[EigenEquation[ϵ],{ϵ,-1}]

Plot[ψ[x][ϵ]/.sol/.eig1,{x,((-L)/2),(L/2)}]
Plot[ψ[x][ϵ]/.sol/.eig2,{x,((-L)/2),(L/2)}]