NDEigensystem potential with a singularity

Consider the following eigenvalue equation

$$-\frac{d^2}{dx^2}\Psi_n(x)+\left(gx^4+\frac{1}{x^2}\right)\Psi_n(x)=E_n\Psi_n(x),\qquad x\in(-10,10),\qquad\Psi(-10)=\Psi(10)=0$$

The boundary of $$x$$ is arbitrary, Ideally it should be $$x\in(-\infty,\infty)$$ but since we are doing numerics we cant impose that. I want to solve this system using NDEigensystem. I have tried the following so far

n = 25;
l = 10;
g = 1/4;
V[x_] := (g*x^4 + 1/x^2);
op = -D[D[\[Psi][x], x], x] + V[x]*\[Psi][x];
bc = DirichletCondition[\[Psi][x] == 0, True];
{vals, funs} = NDEigensystem[{op, bc}, \[Psi][x], {x, -l, l}, n];


When I run this code, I get all sorts of crazy eigenfunctions with singularities. Not ideal at all. But when I run it with a smaller $$l$$, say l=Pi, I get much well behaved functions. How should I adjust my boundaries, if ideally they need to be infinity. Should I somehow consider the potential too?

• When I run this code for l=10 Pi, it runs perfectly fine, and the corresponding solutions look just fine when I plot them, with no singularities. Can you clarify what the issue is, or give an example of when things go wrong? Commented Aug 27, 2021 at 23:37

There are some considerations that should be made before solving the equation.

First of all, your potential, $$V(x)=gx^4+\frac{1}{x^2}$$, is singular at the origin. Hence, an additional condition is needed if you want well-behaved eigenfunctions: vanishing function at $$x=0$$.

Secondly, $$V(x)$$ is even: $$V(x)=V(-x)$$. Therefore, it is sufficient to solve the equation on the semi-infinite interval $$[0,\infty)$$ with boundary conditions $$\Psi_n(0)=\Psi_n(\infty)=0$$. However, in practice, this is not possible as you mention correctly. Thus, you have to truncate your domain and solve on $$[0,l]$$, with $$l$$ sufficiently large.

The natural question: What does it mean sufficiently large $$l$$?

In this particular example, asymptotic analysis provides the answer. It can be easily demonstrated that any eigenfunction at large $$|x|$$ behaves like

$$\Psi_n(x)\sim \exp\left(-\frac{1}{6}|x|^3\ +\ ...\right)\ ,\qquad (g=1/4)$$

see Bender & Orszag for details on a general setting. It means that eigenfunctions decay extremely fast once $$|x|$$ is large. For example, $$l=6$$ leads to $$|\Psi_n(l)|\sim 10^{-16}$$.

Based on your code, I propose the following: replace {x, -l, l} by {x, 0, l} and set l=10. Thus,

n = 25;
l = 10;
g = 1/4;
V[x_] := (g*x^4 + 1/x^2);
op = -D[D[\[Psi][x], x], x] + V[x]*\[Psi][x];
bc = DirichletCondition[\[Psi][x] == 0, True];
{vals, funs} = NDEigensystem[{op, bc}, \[Psi][x], {x,0, l}, n, Method -> {"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];


I inserted some options in Method to increase the accuracy in calculations. Playing a little bit with the value of $$l$$, it seems that the first 10 eigenvalues are reliable when $$l=10$$ and the corresponding eigenfunctions are well-behaved, as expected. See below.

• Can you show the Mathematica code that solves this. This will be much more helpful for the others. Commented May 26, 2023 at 3:36
• As has already been suggested, answers provided by users should be based on Mathematica code. Please take that into consideration and provide the necessary code that answers the question!
– bmf
Commented May 26, 2023 at 7:56
• Thanks for the suggestion @bmf.
– Jam
Commented Jul 20, 2023 at 19:43