Shooting Method for Bound-State Wave Functions

This notebook introduces the shooting method for obtaining wave functions and energies for bound states of the Schrödinger equation.  We consider one-dimensional potentials in the body of the notebook and develop extensions to three dimensions in the exercises that follow.  We begin with a numerical demonstration of how imposition of a boundary condition at infinity leads to a discrete spectrum of energy eigenvalues.  Thus, we write a module for a one-sided shooting method in which we determine the energy that best satisfies a similar boundary condition imposed at a finite radius.  Next we develop a more robust two-sided method which integrates the differential equation inward from both sides and obtains the energy by matching the logarithmic derivative at an interior point.  Our specific example studies the energy splitting between even and odd solutions in a symmetric double-well potential.  The exercises extend these methods to three-dimensional problems.

Initialization

Defaults and packages

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
Off[NDSolve::precw]

$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Utilities`FilterOptions`"];
Needs["Graphics`"]

Introduction

A wide variety of boundary-value problems can be solved using the shooting method in which one formulates a matching condition based upon trial solutions of the differential equation and uses root-finding techniques to adjust a parameter or eigenvalue such that the matching condition is satisfied.  We illustrate this technique by studying the energy eigenvalues for a one-dimensional Schrödinger equation based upon a symmetric double-well potential.  The Schrödinger equation has the form

-((ψ''[x])/(2 m)) + V[x] ψ[x]==e ψ[x]

where m is the mass, V[x] is the potential, e is the energy, and ψ[x] is a wave function which for bound states must satisfy the boundary conditions

(lim_(x -> -∞))ψ[x] == 0
(lim_(x -> +∞))ψ[x] == 
      0

and be everywhere continuous and differentiable.  In general these conditions can be satisfied only for discrete values of the energy, known as eigenvalues or eigenenergies.  Eigenfunctions should be normalized to unit probability according to

(∫_(-∞)^∞)|ψ[
                  x](|^2) dx  = 
      1

For the purposes of illustration, we define the harmonic oscillator potential and a symmetric double-well potential in one dimension that resembles a harmonic oscillator for large x but has a bump in the middle.  The size of the bump has no particular significance.

ho[x_] := ((x^2)/2);
pot[x_] := (x^2) + 5  (E^(-(x^2)))

bottom = FindMinimum[pot[x], {x, 1}];
emin = bottom[[1]]; 
    x0 = x /. bottom[[2]];
{emin, x0}

Plot[pot[x], {x, -5, 5}, PlotRange -> {All, {2, 8}}, 
        Frame -> True, FrameLabel -> {"x", "V[x]"}, 
        PlotLabel -> "a double-well potential"];

Below we use NDSolve to integrate the differential equation with m = 0.5 from xmin=-5 toward  large positive x for a range of energies near one of the eigenvalues.

Clear[f];
f[x_, e_] := 
      Module[{ψ, sol, xmatch = 5}, 
        sol = NDSolve[{-(ψ^″)[x] + 
                    pot[x] ψ[x] == 
                  e ψ[x], ψ[-xmatch] == 
                  0, (ψ^′)[-xmatch] == 
                  (10^(-4))}, ψ, {x, -xmatch, 
                xmatch}][[1]]; 
ψ[x] /. sol]

tmp = Table[f[x, e], {e, 4.43, 4.431, 0.0002}];
Plot[Evaluate[tmp], {x, -5, 5}];

This particular family of solutions possesses one node for moderate x representing approximations to the lowest odd eigenfunction.  For energies near 4.4305 the solutions tend to approach zero for positive x->-xmin, but diverge unless the energy is precisely the correct eigenvalue.   From this figure we deduce that the eigenvalue is between 4.4304 and 4.4306.  A more precise value can be obtained by using FindRoot to locate the energy for which ψ[-xmin]==ψ[xmax]==0.

Note: recent changes to the FindRoot function require the equation to be specified externally with an explicitly numeric argument.

xmatch = 5;
bc[
          e_?NumericQ] := ψ[
            xmatch] /. NDSolve[{-(ψ^″)[x] + 
                    pot[x] ψ[x] == 
                  e ψ[x], ψ[-xmatch] == 
                  0, (ψ^′)[-xmatch] == 
                  (10^(-4))}, ψ, {x, -xmatch, 
                xmatch}][[1]];
evalue = 
      e /. FindRoot[bc[e], {e, 4.43, 4.431}, AccuracyGoal -> 10, 
          WorkingPrecision -> 20]

Despite the 10 digit accuracy we requested for the boundary condition, the solution inevitably diverges if carried too far, either because of the remaining numerical error in the eigenvalue or the accumulated integration error in the solution to the differential equation.

xmin = -5; xmax = 7;
ψ[
          x_] = ψ[
            x] /. NDSolve[{-(ψ^″)[x] + 
                    pot[x] ψ[x] == 
                  evalue ψ[x], ψ[-xmatch] == 
                  0, (ψ^′)[-xmatch] == 
                  (10^(-4))}, ψ, {x, xmin, 
                xmax}][[1]];
Plot[ψ[x], {x, xmin, xmax}];

Clearly one must choose the integration range wisely — if the range is too narrow, inaccuracy in the initial conditions becomes important; if too wide, integration errors accumulate.

Suppose that the potential varies relatively slowly, such that local solutions to the Schrödinger equation can be approximated by

ψ[x] ≈ 
      a (e^(ik x)) + 
        b (e^(-i k x))

where k = (Sqrt[2  m(e - V)]) is the wave number and a and b are arbitrary constants.  In the classically allowed region, where the energy is greater than the local value of the potential, the wave number is real and the convex solutions are oscillatory.  Conversely, in the classically forbidden region, where the energy is smaller than the potential, k = i κ is imaginary and the concave solution is a mixture of exponentially decaying and exponentially growing components.  Eigenfunctions in the classically forbidden region are composed of a delicate balance of growing and decaying solutions whose coefficients vary as the local wave number changes.  In general, numerical integration algorithms encounter trouble when small but inevitable numerical errors upset this delicate balance such that the exponentially growing solution eventually overwhelms the desired eigenfunction.  The classically allowed and forbidden regions for a particular eigenenergy are shown below, along with the corresponding eigenfunction.

tp = x /. 
              FindRoot[
                pot[x] == 
                  evalue, {x, #}] & /@ {-2, - .5,  .5, 
            2};
regions = 
      Graphics[{RGBColor[0.5, 0.05, 0.03], 
          Rectangle[{-5, 
              2}, {tp[[1]], 8}], 
          Rectangle[{tp[[4]], 2}, {5, 
              8}], Rectangle[{tp[[2]], 
              2}, {tp[[3]], 8}], 
          RGBColor[0, 0.26, 0.13], 
          Rectangle[{tp[[1]], 
              2}, {tp[[2]], 8}], 
          Rectangle[{tp[[3]], 
              2}, {tp[[4]], 8}]}]; 
    Plot[{pot[x], evalue, evalue + 10  ψ[x]}, {x, -5, 5}, 
      PlotRange -> {All, {2, 8}}, 
      PlotStyle -> {GrayLevel[0], {GrayLevel[1], Dashing[{0.02, 0.04}]}, 
          RGBColor[0, 0, 1]}, Frame -> True, 
      FrameLabel -> {"x", "V[x], ψ[x]"}, 
      PlotLabel -> "a double-well potential", 
      Prolog -> 
        regions[[1]]];

The classically allowed and forbidden regions are indicated by green or red for the energy shown as a white line.  The corresponding eigenfunction, shifted to the energy level and with arbitrary normalization, is shown in blue.

A more robust version of the shooting method attempts to match

(Subscript[ψ, 1])[(Subscript[x, m])] == (Subscript[ψ, 2])[(Subscript[x, m])]
(Subscript[ψ, 1]^′)[(Subscript[x, m])] == (Subscript[ψ, 2]^′)[
        (Subscript[x, m])]

where (Subscript[ψ, 1]) is obtained by integration from xmin to xmatch while (Subscript[ψ, 2]) is integrated from xmax to xmatch where xmin<xmatch<xmax.  Equivalently, if we match the two values of the logarithmic derivative (ψ^′)[x]/ψ[x] at an appropriately chosen interior match point, it is not necessary to normalize the solution during the search for the eigenenergy.

Unfortunately, it is difficult or impossible to formulate rules that can be used to choose the starting conditions and search intervals in a general automatic procedure that will work for any potential.  Thus, it is necessary to tailor the method of solution to the special characteristics of the problem at hand, which requires the kind of insight that only experience can provide.  In this notebook we examine in some detail the solutions to a particular double-well potential that is sufficiently challenging to illustrate some of the issues that must be faced yet simple enough to solve without more specialized techniques.

Implementation

One-sided shooting method

The following function solves the eigenvalue problem by integrating the Schrödinger equation from xmin to xmax and finding an energy for which ψ[xmax]==0.

ShootBSWF1[pot_, m_, emin_, emax_, xmin_, xmax_, xmatch_, 
        opts___Rule] := 
Module[{ψ, x, y, trial, bc, 
          evalue, e, norm, optsRoot, optsND, 
          opts2}, 
optsRoot = 
          FilterOptions[FindRoot, opts]; 
        optsND = 
          FilterOptions[NDSolve, opts]; 
trial[
            e_] := NDSolve[{-(y^″)[x]/(2  m) + 
                    pot[x] y[x] == e y[x], 
                y[-xmatch] == 
                  0, (y^′)[-xmatch] == (10^(-6))}, 
              y, {x, xmin, xmax}, 
              optsND][[1]]; 
bc[e_?NumericQ] := 
          y[xmatch] /. trial[e]; 
evalue = 
          If[emin == emax, emin, 
            e /. With[{opts2 = optsRoot}, 
                FindRoot[bc[e] == 0, {e, emin, emax}, 
                  opts2]]]; 
ψ = 
          y /. trial[evalue]; 
norm = 
          1/(Sqrt[NIntegrate])[ψ[x]^2, {x, xmin, xmax}]; 
{energy -> evalue, 
          wavefunction -> (norm ψ[#] &), 
          normalization -> norm}]

solutions = 
        wavefunction[x] /. 
          Table[ShootBSWF1[ho, 1, n + 0.4, n + 0.6, -5, 5, 5, 
              MaxIterations -> 50], {n, 0, 3}];

Plot[Evaluate[solutions], {x, -5, 5}];

Interior matching condition

The following function compares values of the logarithmic derivative (ψ^′)[x]/ψ[x] at an interior point xmatch obtained by integrating either from xmin to xmatch or from xmax to xmatch.  The difference, or mismatch, between these values can then be used to determine the eigenenergy.

mismatch[e_, m_, pot_, xmin_, xmax_, xmatch_, 
        opts___Rule] := 
Module[{bc1, bc2, x, 
          y}, 
bc1 = (y^′)[xmatch]/
              y[xmatch] /. (NDSolve[{-(y^″)[x]/(2  
                                m) + pot[x] y[x] == e y[x], 
                    y[xmin] == 0, (y^′)[xmin] == 
                      (10^(-6))}, y, {x, xmin, xmatch}, 
                  opts][[1]]); 
bc2 = (y^′)[xmatch]/
              y[xmatch] /. (NDSolve[{-(y^″)[x]/(2  
                                m) + pot[x] y[x] == e y[x], 
                    y[xmax] == 0, (y^′)[xmax] == 
                      (10^(-6))}, y, {x, xmax, xmatch}, 
                  opts][[1]]); 
        bc1 - bc2]

Double-sided shooting method

The following function integrates from xmin to xmatch and then from xmax to xmatch and uses FindRoot to determine an energy for which the mismatch between two logarithmic derivative vanishes.  The two sides of the solution are then normalized and patched together to generate the eigenfunction.

ShootBSWF2[pot_, m_, emin_, emax_, xmin_, xmax_, xmatch_, 
        opts___Rule] := 
Module[{ψ, x, y, x0, trial, 
          eq, bc, bc1, bc2, ψ1, ψ2, ratio, evalue, e, norm1, norm2, 
          norm, optsRoot, optsND, opts2}, 
optsRoot = 
          FilterOptions[FindRoot, opts]; 
        optsND = 
          FilterOptions[NDSolve, opts]; 
eq[e_, x_, 
            x0_] := {-(y^″)[x]/(2  m) + 
                pot[x] y[x] == e y[x], 
            y[x0] == 0, (y^′)[x0] == 
              (10^(-6))}; 
y1[e_, x_] := 
          y[x] /. NDSolve[eq[e, x, xmin], y, {x, xmin, xmatch}, 
                optsND][[1]]; 
bc1[e_] := 
          D[y1[e, x], x]/y1[e, x] /. x -> xmatch; 
y2[
            e_, x_] := 
          y[x] /. NDSolve[eq[e, x, xmax], y, {x, xmax, xmatch}, 
                optsND][[1]]; 
bc2[e_] := 
          D[y2[e, x], x]/y2[e, x] /. x -> xmatch; 
bc[
            e_?NumericQ] := bc1[e] - bc2[e]; 
evalue = 
          If[emax == emin, emin, 
            e /. With[{opts2 = optsRoot}, 
                FindRoot[bc[e], {e, emin, emax}, 
                  opts2]]]; 
ψ1[x_] = 
          y1[evalue, x];  ψ2[x_] = 
          y2[evalue, 
            x]; 
ratio = ψ1[xmatch]/ψ2[
              xmatch]; 
norm1 = 
          NIntegrate[ψ1[x]^2, {x, xmin, 
              xmatch}]; 
norm2 = 
          NIntegrate[ψ2[x]^2, {x, xmatch, 
              xmax}]; 
norm = 
          1/(Sqrt[norm1  +  (ratio^2) norm2]); 
ψ[
            x_] = If[x ≤ xmatch, norm ψ1[x], 
            norm ratio ψ2[x]]; 
{energy -> 
            evalue, wavefunction -> (ψ[#] &), 
          normalization -> norm}]

A plotting function

The following function plots the potential, eigenenergy, and eigenfunction together.

Clear[PlotBSWF];
PlotBSWF[solution_, xmin_, xmax_, opts___Rule] := 
        Module[{e = energy /. solution, ψ, 
            f}, 
ψ[x_] = 
            wavefunction[x] /. solution; 
f[x_] = 
            e + If[xmin < x < xmax, ψ[x], 
                0]; 
Plot[{pot[x], f[x], e}, {x, xmin, 
              xmax}, opts, 
            PlotStyle -> {RGBColor[1, 0, 0], Dashing[{0}], 
                Dashing[{0.02, 0.04}]}, Frame -> True, 
            FrameLabel -> {"x", "V, ψ"}, 
            PlotLabel -> "Wave Function with Energy = " <> 
                ToString[N[energy /. solution, 4]]]];

Test case: harmonic oscillator

The eigenvalue problem for the harmonic oscillator potential V[x] = (1/2) k (x^2) can be solved symbolically, with exact eigenenergies (Subscript[ε, n]) = (n + (1/2)) ω where ω = (Sqrt[k/m]) is the classical oscillation frequency.  Using the double-sided shooting method, the energy eigenvalues are determined by matching the logarithmic derivative of the wave functions obtained by integration from either side of the potential toward the interior.  Although the match point can be chosen almost anywhere, for a symmetric potential it is natural to choose the center of the potential.  Even solutions have finite wave functions but vanishing derivatives, so that the mismatch between logarithmic derivatives exhibits a simple zero.  Odd solutions have finite derivatives but vanishing wave functions, so that the mismatch exhibits a singularity at the eigenvalue.The figure below illustrates the matching condition used by the symmetric shooting method.  Note that we must increase the default AccuracyGoal to obtain a smooth curve for higher energies.

Plot[
        mismatch[e, 1, ho, -5, 5, 0, AccuracyGoal -> 10], {e, 0, 
          10}, PlotRange -> {All, {-10, 10}}, 
        PlotPoints -> 100, Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=1"];

Although we matched at the center above, that point is not suitable for FindRoot because the singularity for odd solutions makes the secant method unstable.  For construction of the wave functions it is better to use a small value of x≠ 0 where a wave function node is unlikely.  The lowest four eigenfunctions are displayed below using this method and are compared with exact wavefunctions.  Note that the exact wave functions employ the parameter α = (Sqrt[m k]) and differ in phase from the numerical wave functions by (((-))^n).

ψ[n_, α_] = 
        Function[
          x, (Sqrt[(α/((2^n) n! (Sqrt[π])))]) 
            Exp[-(α^2) (x^2)/2] 
            HermiteH[n, α x]];

solutions = 
        wavefunction[x] /. 
          Table[ShootBSWF2[ho, 1, n + 0.4, n + 0.6, -5, 5, 0.5, 
              MaxIterations -> 50], {n, 0, 3}];

\(\(DisplayTogether[
        Plot[Evaluate[solutions], {x, \(-5\), 5}, Frame -> True, 
          FrameLabel -> {\*""\!\(\@\(m\\ k\)\)x"", \*""(Subscript[ψ, n])""}, PlotLabel -> "points: exact, curves: numerical"], 
        ListPlot[
          Flatten[Table[{x, \(\((\(-1\))\)\^n\) \(ψ[n, 1]\)[
                    x]}, {x, \(-5\), 5, 0.1}, {n, 0, 3}], 1], 
          PlotStyle -> AbsolutePointSize[2]]];\)\)

Properties of the double-well potential

For low energies, the potential plotted below can be described as two wells separated by a barrier.

FilledPlot[{pot[x], 4}, {x, -5, 5}, 
        PlotRange -> {All, {2, 8}}, Frame -> True, 
        FrameLabel -> {"x", "V[x]"}, 
        PlotLabel -> "a double-well potential"];

Suppose that the energy is between the bottom of the well and its central maximum, 2.6 < ε < 5.  There will then be two pairs of classical turning points, symmetrically placed, where the energy and potential are equal, such that ε == V[(Subscript[x, tp])].  In the outer regions and the central region where the potential is above the energy, the wave function is a mixture of exponentially decaying and exponentially growing components, whereas it becomes oscillatory within the two wells in the regions where the energy is below the potential.   If the energy is near the top of the barrier, the wave function is not strongly suppressed within the barrier, but if the energy is very close to the bottom of the well, the central barrier then looks very thick and there is very little penetration of the wave function into this classically forbidden region.  We will find that the thick barrier case is numerically very challenging.  The effective thickness of the barrier is determined by the mass.

It is useful to construct an equivalent harmonic oscillator that is a good approximation near the bottom of a well.

potseries = Series[pot[x], {x, x0, 2}]

pot1[x_] = Normal[potseries]

Plot[{pot[x], pot1[x]}, {x, -3, 3}, 
        PlotStyle -> {Dashing[{}], Dashing[{0.02, 0.04}]}];

(Subscript[k, eff]) = 2  SeriesCoefficient[potseries, 2]

ω = (Sqrt[(Subscript[k, eff])/m])

For small ω we can approximate low-lying energy levels as follows.

(Subscript[ε, n])_ = emin + (n + (1/2)) ω

Thus, roughly speaking, the thin barrier applies when m ≤ 1 and the thick barrier when m>10.

Narrow barrier

Low-lying states

Approximating each well by an independent harmonic oscillator, the lowest energy levels for m = (1/2) would be

Table[(Subscript[ε, n]) /. m -> (1/2), {n, 0, 1}]

Because the wave functions for each well overlap within the barrier, the even and odd combinations of these wave functions have different energies, with the even combination being lower because the kinetic energy contribution is smaller.  Thus, the single-well eigenvalues separate into pairs and the splitting between the even and odd partners increases with n because the higher wave functions overlap more strongly in the central region.  The figure below reveals even solutions at e≈3.89 and 6.55 and odd solutions at e ≈ 4.43 and 8.13.

Plot[mismatch[e, 0.5, pot, -5, 5, 0], {e, 3.5, 10}, 
        PlotRange -> {All, {-10, 10}}, PlotPoints -> 100, 
        Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=0.5"];

However, matching at the center does not work particularly well for construction of the wave functions, especially for odd solutions.  Furthermore, it is also difficult to ensure that the FindRoot will settle upon the desired solution.  Therefore, a certain amount of trial and error is needed with these algorithms to choose parameters which will produce the desired solutions.  The following choices produce the lowest four solutions for this case.

solutions = {
ShootBSWF2[pot, 0.5, 3, 4, -5, 5, 
            0.5, MaxIterations -> 50], 
          ShootBSWF2[pot, 0.5, 4.4, 4.5, -5, 5, 0.5, 
            MaxIterations -> 50], 
ShootBSWF2[pot, 
            0.5, 6, 7, -5, 5, 0.5, 
            MaxIterations -> 50], 
ShootBSWF2[pot, 
            0.5, 8, 9, -5, 5, 0.5, 
            MaxIterations -> 50]};

plot1 = 
        PlotBSWF[solutions[[1]], -5, 
          5, PlotRange -> {All, {2.5, 5.5}}, 
          DisplayFunction -> Identity];
plot2 = 
        PlotBSWF[solutions[[2]], -5, 
          5, PlotRange -> {All, {2.5, 5.5}}, 
          DisplayFunction -> Identity];
plot3 = 
      PlotBSWF[solutions[[3]], -5, 
        5, PlotRange -> {All, {5.5, 7.5}}, 
        DisplayFunction -> Identity]; 
    plot4 = PlotBSWF[
        solutions[[4]], -5, 5, 
        PlotRange -> {All, {7.5, 9}}, 
        DisplayFunction -> Identity];
plots = {{plot1, plot2}, {plot3, plot4}};

Show[GraphicsArray[plots], 
        DisplayFunction -> $DisplayFunction];

Plot[
        Evaluate[wavefunction[x] /. solutions], {x, -5, 5}];

Higher energies

Higher energies must be integrated over larger ranges to ensure that the initial conditions are adequate.  As the energy increases, the central bump becomes less important and the spacing between even and odd levels approaches the uniform spacing of the simple harmonic oscillator, but now with k = 1.

Plot[
        mismatch[e, 0.5, pot, -7, 7, 0, MaxSteps -> 5000], {e, 3, 
          22}, PlotRange -> {All, {-50, 50}}, 
        PlotPoints -> 100, Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=0.5"];

sol = 
        ShootBSWF2[pot, 0.5, 21, 22, -7, 7, 0.5, 
          MaxIterations -> 50];
PlotBSWF[sol, -7, 7, 
        PlotRange -> {All, {21, 22.3}}];

Thick barrier

Next, consider a thick-barrier case with m = 20.  It is useful to estimate the ground-state energy using the harmonic oscillator approximation for an isolated well.  This energy is quite close to the bottom of well, so that the barrier appears to be rather thick.  Nevertheless, we might expect the true energy to be slightly smaller than this estimate because our model potential rises higher than the true potential on one side.

(Subscript[ε, 0]) /. m -> 20

Although there is a hint of glitch near 2.88, the zero/pole structure is so sharp that it is easily missed unless a very large number of points is used by the plotting routine.

Plot[
        mismatch[e, 20, pot, -3, 3, 0, MaxSteps -> 5000], {e, emin, 
          3}, PlotRange -> {All, {-50, 50}}, PlotPoints -> 500, 
        Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=20"];

One can zoom in upon the glitch, but it is difficult to find because the even and odd solutions are so close together.

Plot[
        mismatch[e, 20, pot, -3, 3, 0, MaxSteps -> 5000], {e, 
          2.88715, 2.8872}, PlotRange -> {All, {-50, 50}}, 
        PlotPoints -> 500, Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=20"];

Recognizing that the even and odd solutions are practically degenerate, one can instead locate the energies using just one side of the potential and choosing a match point at any convenient point.  However, in this case only the zeros of the matching condition represent true solutions — unless we accidentally select a match point that coincides precisely with a node in the wave function, the poles are spurious.  By starting the wave function at zero, this method finds the odd solutions only, but we know that for each odd solution there is an even solution at slightly smaller energy.   Note we need to increase AccuracyGoal to obtain a clean plot.

Plot[
        mismatch[e, 20, pot, 0, 3, 2, AccuracyGoal -> 10, 
          MaxSteps -> 5000], {e, 2.85, 5}, 
        PlotRange -> {All, {-50, 50}}, PlotPoints -> 100, 
        Frame -> True, 
        FrameLabel -> {"energy", "mismatch"}, 
        PlotLabel -> "matching condition for m=20"];

The eigenvalues for the odd solutions are now readily apparent.  We locate the lowest odd solution using ShootBSWF1 for one side of the potential.

sol1 = 
      ShootBSWF1[pot, 20, 2.8871, 2.88715, 0, 3, 3, AccuracyGoal -> 10, 
        WorkingPrecision -> 32, MaxSteps -> 10000, 
        MaxIterations -> 200]

PlotBSWF[sol1, 0, 3];

Then searching for a slightly lower eigenvalue, we obtain the ground-state solution below.

sol = 
      ShootBSWF2[pot, 20, 2.8871, 2.88715, -3, 3, 0, 
        AccuracyGoal -> 10, WorkingPrecision -> 32, 
        MaxIterations -> 200]

PlotBSWF[sol, -3, 3];

The splitting between the even and odd solutions is extremely small and becomes even smaller as the mass increases.

(energy /. sol1)  - (energy /. sol)

Clearly, it becomes rather difficult with this method to obtain an accurate value for the ground-state splitting as the mass becomes large.

\(TraditionalForm\`\(-\(∇\^2 ψ[\(r\& ⇀ \)]\/\(2\ m\)\)\)\  + \ V[r]\ ψ[\(r\& ⇀ \)] == 
      e\ ψ[\(r\& ⇀ \)]\)

\(TraditionalForm\`ψ[r\&⇀] == \ 
      R\_\(n, ℓ\)[r] Y\_\(ℓ, m\)[r\&^]\)

(Subscript[R, n, ℓ]^″)[r] + (2/r) 
          (Subscript[R, n, ℓ]^′)[
            r]  - ((ℓ (ℓ + 
                    1))/( (r^2))) 
          (Subscript[R, n, ℓ])[
            r] + ((2  m)/(ℏ^2)) ((Subscript[E, n, ℓ]) - 
              V[r]) (Subscript[R, n, ℓ])[r] == 0

\(TraditionalForm\`r -> 
      0\ \ ⟹\ \ R\_\(n, ℓ\)\  ∝ \ r\^ℓ\)

\(TraditionalForm\`ℓ = \(0\ \ ⟹\ \ R\_\(n, ℓ\)[r\_0] = 1\), \ \ R\_\(n, ℓ\)\%′[
        0] = \(-ε\)\)
ℓ > 
        0  ⟹  (Subscript[R, n, ℓ])[(Subscript[r, 0])] = 
      0,   (Subscript[R, n, ℓ]^′)[
        0] = +ε

(∫_0^∞) (Subscript[R, n, ℓ])[
                r]^2 (r^2) dr  =  
      1

(Subscript[ψ, n, ℓ])[
        r]  =  (Sqrt[((2  Γ[
                    n - ℓ + 1])/Γ)[n + (3/2)]]) 
        Exp[-m ω (r^2)/2] (r^ℓ) 
        (Subscript[L, n - ℓ]^((2  ℓ + 1)/2))[
          m ω (r^2)] ,     (Subscript[E, n, ℓ]) = (2  
              n - ℓ + (3/2)) ℏ ω

\(TraditionalForm\`ψ\_\(n, ℓ\)[
        r]\  = \ \(\(2  κ\^\(3/2\)\)\/n\^2\) \(\@\(\(\((n - ℓ + 1)\)!\)\/\(\((n + ℓ)\)!\)\)\) 
        Exp[\(-\(\(κ\ r\)\/n\)\)] \(\((\(2  κ\ r\)\/n)\)\^ℓ\) 
        L\_\(n - ℓ - 1\)\%\(2  ℓ + 1\)[\(2  κ\ r\)\/n]\ , \ \ \ \ E\_\(n, ℓ\) = \(-\(E\_0\/n\^2\)\)\)

\(TraditionalForm\`V[r]\  = \ 
      V\_C\/\(Exp[\((r - r\_C)\)/a\_C] + 1\)\ \  + \ \(ℏ\^2\/\(\(m\_π\^2\) 
                c\^2\)\) \(1\/r\) \(∂\/∂r\) \(V\_LS\/\(Exp[\((r - r\_LS)\)/a\_LS] + 
                1\)\) \(L\& ⇀ \)·\(σ\& ⇀ \)\)

\(⟨\(L\& ⇀ \)·\(σ\& ⇀ \)⟩\) = {({{ℓ, j = ℓ + (1/2)}, {-ℓ - 1, j = ℓ - (1/2)}})

V[r]  = (Subscript[V, 0]) ( ((((r/(Subscript[r, 0]))))^4) - 
            2 ((((r/(Subscript[r, 0]))))^2) )