代码是
Clear["Global`*"]
m := 1.320;
\[Sigma] := 0.200;
\[Alpha]s := 0.35;
\[Mu]R := m/4;
\[Delta] := 10^-10;
V[r_] := -4/3*\[Alpha]s/r + \[Sigma]*r;
Veff[r_] := 2*\[Mu]R*V[r];
SolveRadialEq[\[CurlyEpsilon]_] :=
NDSolve[{u''[r] == (Veff[r] - \[CurlyEpsilon]^2/4 + m^2) *u[r],
u[0] == \[Delta], u'[0] == 1}, u, {r, 0, 20}];
CheckBoundaryCondition[sol_, \[CurlyEpsilon]_] :=
Module[{uSol, value}, uSol := u[r] /. First[sol]; value := uSol[20];
value];
FindEnergyLevel[l_, \[CurlyEpsilon]Min_, \[CurlyEpsilon]Max_,
tolerance_] :=
Module[{\[CurlyEpsilon]Mid, sol, check},
While[\[CurlyEpsilon]Max - \[CurlyEpsilon]Min >
tolerance, \[CurlyEpsilon]Mid := (\[CurlyEpsilon]Min + \
\[CurlyEpsilon]Max)/2; sol := SolveRadialEq[l, \[CurlyEpsilon]Mid];
check := CheckBoundaryCondition[sol, \[CurlyEpsilon]Mid];
If[check >
0, \[CurlyEpsilon]Max = \[CurlyEpsilon]Mid, \[CurlyEpsilon]Min \
= \[CurlyEpsilon]Mid];
];
\[CurlyEpsilon]Mid];
\[CurlyEpsilon]Ground = FindEnergyLevel[0, -1, 1, 10^-6]
