I tried solving the eigenvalue problem of a 2nd-order ODE $$[b^2(k-2)^2y^2-2b(k-2)(1+2ky)+4k^2+b^2(k-2)3y]f(y) \\- 3b(3by-2)f'(y)\\-(3by-2)^2f''(y)=\lambda f(y)$$ with NDEigenvalue. $b,k$ are just real parameters. Note that the long bracket in front of $f(y)$ is just a quadratic polynomial of $y$. Let's say $b>0$ for simplicity. The boundary condition is Dirichlet at infinity or at least very far away to get a completely decayed convergent solution and here I just use Dirichlet at empirically large enough yL and yR.
Thanks to @bbgodfrey, who reminded me of a singularity from $(3by_\mathrm{sgl}-2)=0$. I updated much of this post and prefer not showing the very naive old version.
- It appears from the equation that the singularity renders both the 1st- and 2nd-order derivatives to contribute less and less when approaching ${y_\mathrm{sgl}}-$, beyond which the solution (possibly??) varies little. Numerically, this holds nice when $b$ is small while the turning point will be smoothed and shifted further rightward when $b$ is larger (but why?).
- Numerically, another important and intriguing finding is that increasing $b$ has the effect of squeezing the solution towards $y_\mathrm{sgl}$, which is not yet intuitively clear to me. Therefore, the larger $b$ becomes, the more unstable the numerical solution might be, I guess.
- But what is the solution to be squeezed? Actually, a parent and simpler equation is $$[b^2(k-2)^2y^2-2b(k-2)(1+2ky)+4k^2]f(y) \\-4f''(y)=\lambda f(y),$$ whose $n$th solution is basically related to the $n$th Hermite polynomial with $n-1$ node(s). From the figures below, one would tentatively guess that it's this to be squeezed. Rather surprising to me, the node structure looks intact.
When $b$ is small like 0.03, the solution seems probably not that bad. But for larger $b$ like 0.5,1.0 or so, the resultant eigenvalues are very unstable against tuning the MeshOptions. Unfortunately, I don't know if there is an analytical solution. So I hope there could be some way around or even other methods to improve or reassure and it certainly also helps if one can somewhat understand the numerical observations.
k = 0; b = 0.03; Nless = 20;
ysgl = Solve[(-2 + 3 b y) == 0, y][[1, 1]] // Last; yR = 1.3 ysgl; yL = -Sign[b] 7/Sqrt[Abs@b];
lhs = (b^2 (-2 + k) y (3 + (-2 + k) y) - 2 b (-2 + k) (1 + 2 k y) + 4 k^2) f[y]
- (-2 + 3 b y) (3 b D[f[y], y] + (-2 + 3 b y) D[f[y], {y, 2}]);
bc = DirichletCondition[f[y] == 0, True];
{vals, funcs} =
NDEigensystem[{lhs, bc}, f[y], {y, yL, yR}, Nless,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
-> {"MaxCellMeasure" -> 0.01, "MeshOrder" -> 2}}}}];
vals
The Norm of first ten solutions plotted when $b=0.03$.
when $b=0.06$
Update of the analytic solution and a puzzling point
It's really nice to have an analytic solution like the one by @bbgodfrey. We can actually first shift $y\rightarrow y+\frac{2}{3b}$ in the ODE. Then the singularity $y_\mathrm{sgl}=0$. The solution is $f(y)=\mathrm{e}^{-\frac{1}{3}(k-2)y}y^{\frac{\Delta}{9b}}u(y)$ where $\Delta=\sqrt{16(k+1)^2-9\lambda}$ and $u$ is polynomial.
However, the exponential factor in $f(y)$ will diverge on the right of $y_\mathrm{sgl}$. It looks that analytic solution doesn't allow a physical solution vanishing at $+\infty$, although it's good if one sets Dirichlet b.c. exactly at $y_\mathrm{sgl}$. But I would expect a solution even for $y>y_\mathrm{sgl}$. If not, this is really surprising since the original physical meaning of y is merely a coordinate in $(−\infty,+\infty)$, not anything like a radius defined in $(0,+\infty)$.
On the other hand, the numerical solution always looks vanishing for $y>y_\mathrm{sgl}$. From the analytic solution, @bbgodfrey claims an upper bound of eigenvalues and higher ones are spurious. What I observe in numerics is twofold. For those 'good' eigenvalues, $f(y_\mathrm{sgl})=0$ seems sharp. For those 'spurious' eigenvalues, $f(y_\mathrm{sgl})=0$ or a turning point seems smoothed. I really have no idea if numeric solution for $y>y_\mathrm{sgl}$ is just an artifact or anything else.
(-2 + 3 b y) == 0, and this is the source of your problem, I believe., $\endgroup$