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.

  • $\begingroup$ What version of Mathematica do you use? I use version 12 for macos, and running the code from your first code block works just fine... $\endgroup$ Commented Dec 23, 2020 at 15:09
  • $\begingroup$ I have 11.2. So you run it, without warnings and getting the $-3$ eigenvalue as well? $\endgroup$
    – Smerdjakov
    Commented Dec 23, 2020 at 15:17
  • $\begingroup$ Yes, it runs fine without any warnings or errors. And the smallest numerical eigenvalue seems to be -2.90075 which is not very far from -3. (Certainly this will come closer to -3 if you refine the discretization.) $\endgroup$ Commented Dec 23, 2020 at 15:20
  • $\begingroup$ Thanks for this. The question is then, is there anything that could be done on 11.2 to get there? This is puzzling though, I wonder if the same underlying reason results in the second method I used failing. $\endgroup$
    – Smerdjakov
    Commented Dec 23, 2020 at 15:24
  • 1
    $\begingroup$ @Smerdjakov I am also got eigenvalue of -2.97433 with v.12.0.0 and without bc ` DirichletCondition[w[y] == 0, True]` at L=10. So let try L=10 and remove bc and option Method -> {"Arnoldi", "Criteria" -> "RealPart"} in v. 12.2 $\endgroup$ Commented Dec 23, 2020 at 15:39

1 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):

    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Load the package:


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}]

Plot of the Evans function

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]

enter image description here

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.

  • 1
    $\begingroup$ thanks for this. I found a way too to solve these equations, but I will certainly look at your tool for a validation if not anything else, thanks a lot. I will check it tomorrow first thing. $\endgroup$
    – Smerdjakov
    Commented Dec 29, 2020 at 21:23
  • $\begingroup$ thanks for pointing out this nice package, it works very well, much more elegant than my caveman approach. I am going through the theory now, the paper you suggested on your GitHub page. Just a question, is the package expected to work if the potential $U$ in the definition is not given in closed-form, but as an interpolating function coming out of a NDSolveValue solution? I get errors but I am not so sure if they are package-related or else. The code in your repository is out of my current Mathematica understanding so I cannot check it. $\endgroup$
    – Smerdjakov
    Commented Dec 30, 2020 at 8:36
  • 1
    $\begingroup$ Yeah, my code is somewhat archaic. In theory it should work with an interpolation function, although it looks like that breaks the decaying boundary conditions at the moment, because I use Limit. That is fixable though. If you try to put in the zero derivative at 0 and zero value at L, then it should work with the interpolation function $\endgroup$
    – SPPearce
    Commented Dec 30, 2020 at 13:36
  • $\begingroup$ that is great, your example notebook is also very helpful. Thanks a lot $\endgroup$
    – Smerdjakov
    Commented Dec 30, 2020 at 17:58
  • 1
    $\begingroup$ You're welcome, I hope it is useful to you, feel free to continue to ask me questions about it. I have not got round to finishing off a paper on the method, if only I had the time... $\endgroup$
    – SPPearce
    Commented Dec 30, 2020 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.