Numerical solution of Eigenvalue equation

Question

I am interested in finding eigenvalues of Schrödinger-type equations, a prototype example being

$$- w^{\prime \prime } (y) - 6 \operatorname{sech}^2 (y) w(y) + w(y) = \lambda w(y)$$

I posted a question on the same equation in Solving a Sturm-Liouville problem and got great help there, both on symbolic methods as well as a numeric approach which delivers the sought result.

Out of curiosity, I have experimented also NDEigensystem and another alternative, described below, and would like to know what is that I am doing wrong with those, as I am not able to recover the same results.

Using NDEigensystem, setting Dirichlet boundary conditions at $\pm 20$ to mimick the real conditions at infinity, $w(y) \to 0$ as $ y \to \infty$ I tried

 L = 20
{eigN2, funcN2} = 
N @ DEigensystem[{w[y] - 6 Sech[y]^2 w[y] - w''[y], 
 DirichletCondition[w[y] == 0, True]}, w[y], {y, -L, L}, 12, 
 Method -> {"Arnoldi", "Criteria" -> "RealPart"}];
 Plot[funcN2, {y, -20, 20}, PlotLegends -> eigN2]

I get a warning message

      ParametricNDSolveValue::bvluc: The equations derived from the boundary 
        conditions are numerically ill-conditioned. The boundary conditions may 
      not be sufficient to uniquely define a solution. If a solution is computed, 
        it may match the boundary conditions poorly. 

Is this a working precision issue? I found a post stating that NDEigensystem cannot work with arbitrary precision. Can the issue be circumvented by a better choice of boundary conditions or what else?

In any case, I get only positive eigenvalues (which do exist and are related to oscillatory eigenfunctions, is the boundary conditions definition not satisfactory?). I used

Method -> {"Arnoldi", "Criteria" -> "RealPart"}

as I think it should sort the eigenvalues by magnitude, as I found on another post, but I am not successful (there should be an eigenvalue equal to $-3$), as documented on my first post on the equation linked above. The equation as written here is obtained by a change of variables from the original one, as documented in the answer by bbgodfrey.

I have also tried the approach detailed in the following post Find eigen energies of time-independent Schrödinger equation. I know the eigenfunction I am interested is even, so I use the b.c. $w'(0) = 0$:

U[x_] := -6*Sech[x]^2
L = 20
pfun2 = ParametricNDSolveValue[{ -ψ''[x] + 
 U[x] ψ[x] + ψ[x] == Ei ψ[x], ψ'[0.] == 
0., ψ[L] == 0.}, ψ, {x, -L, L}, {Ei}]
Show[Plot[U[x], {x, -5, 5}, PlotStyle -> Directive[Thick, Green]], 
Plot[Evaluate@Table[Re[pfun2[Ei][x]], {Ei, -3, 2, 1}], {x, -L, L}], 
 PlotRange -> {All, {-6, 10}}]
val2 = Map[FindRoot[pfun2[Ei][-L], {Ei, #}] &, {1, 2, 3}]
   {{Ei -> 1.02897}, {Ei -> 2.00241}, {Ei -> 3.00088}}

and there is no negative eigenvalue, nor the related eigenfunction.

Any hint would be greatly appreciated.

Answer 1

Ah, this kind of $\operatorname{sech}^2$-based equations with decaying boundary conditions are familiar to me, as they cropped up during during my PhD (particularly this article),

And I happen to have written a package to solve eigenvalue equations of this type, using the method that I used during my PhD. This uses the Evans function, and these kind of conditions are included in my implementation.

The Evans function is an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, my github page and this PDF, or search for CompoundMatrixMethod to see my previous answers here.

Install the package (also available on my github page):

Needs["PacletManager`"]
    PacletInstall["CompoundMatrixMethod", 
    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Load the package:

Needs["CompoundMatrixMethod`"] 
Clear[L]

First we need to convert the system into the appropriate form:

U[x_] := -6*Sech[x]^2
sys[L_] = ToMatrixSystem[-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], {}, ψ, {x, -L, L}, Ei]

where I haven't specified a domain length or any boundary conditions; the way I have the package set up is to assume that the solution should decay if no boundary conditions are given (see the notebook of examples in the github).

Now we can evaluate the Evans function for a given value of $Ei$ (and domain length $L$), say $Ei=-5$ and $L=10$:

Evans[-5, sys[10]]
(* 0.424581 *)

This is not zero, so $Ei=-5$ is not an eigenvalue of the system. So we can use FindRoot:

FindRoot[Evans[Ei, sys[10]], {Ei, -5}]
{Ei -> -3.}

If we plot the Evans function, we can see that it is smooth and has the root at $Ei=-3$.

Plot[Evans[Ei, sys[10]], {Ei, -5, 0}]

Note that the values of this function are around 1, as the Compound Matrix Method that my package uses helps to remove the exponential growth that otherwise is present.

If you vary L, you can see that the eigenvalue converges to -3 when you get to about $L=5$:

Ei /. Table[FindRoot[Evans[Ei, sys[L]], {Ei, -5}], {L, 1, 10}]
(* {-3.40324, -3.01715, -3.00035, -3.00001, -3., -3., -3., -3., -3., -3.} *) 

This doesn't give you the eigenvectors immediately, I have not worked out an easy way to code that into the package. But for this second order equation you can integrate from zero to L with your one boundary condition (that it is even) and an arbitrary value:

eigfun = NDSolveValue[{-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], ψ[0] == 1, ψ'[0] == 0} /.  Ei -> (-3), ψ, {x, 0, 10}]
Plot[eigfun[x], {x, 0, 5}, PlotRange -> All]

Regarding including numerical interpolation functions, my code with the decaying BCs doesn't currently work for non-analytic functions, but if you give boundary conditions it works fine:

uu = FunctionInterpolation[U[x], {x, -10, 10}, 
    InterpolationPoints -> 200, InterpolationOrder -> 5] // Quiet;
sysExact[L_] = ToMatrixSystem[-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], {ψ'[0] == 0, ψ[L] == 0}, ψ, {x, 0, L}, Ei];
sysIntpol[L_] = ToMatrixSystem[-ψ''[x] + uu[x] ψ[x] + ψ[x] == Ei ψ[x], {ψ'[0] == 0, ψ[L] == 0}, ψ, {x, 0, L}, Ei];

These give the same result:

FindRoot[Evans[Ei, sysExact[10]], {Ei, -5}]
FindRoot[Evans[Ei, sysIntpol[10]], {Ei, -5}]

If you have any further questions please just ask.