Consider an ODE $$ \psi''(x) + \left[a^2 - x^{2n}\right]\psi(x) = 0 $$ for $n\in \mathbb{Z}$ and $a\in \mathbb{R}$. The problem originally came from QM and is about finding energy levels in a given potential, so I'm looking for regular solutions (and specific $a$'s which correspond to them) which vanish at $x = \pm \infty$.
So far I've attempted a substitution $\psi = \exp\left(-\frac{x^{n+1}}{n+1}\right) u(x)$ which leads to
$$ u''(x) - 2x^n u'(x) + [a^2-nx^{n-1}] u(x)=0 $$
then I used $x^{n+1} = t$ and got to
$$ (n+1)^2t^2u'' +\bigl[(n+1)nt-2(n+1)t^2\bigr]u' + \bigl[a^2 t^{\frac{2}{n+1}} - n t\bigr]u = 0 $$
Which is "almost nice" apart from the non-integer power in the $a^2$ term. I do not have an idea how can I proceed further.
Another approach I've tried factorizing the initial differential operator
$$ \hat{H} = -\frac{d^2}{dx^2} + x^{2n} = A^*A + \epsilon $$
with
$$ A = \frac{d}{dx} + \beta(x),\;\; A^* = -\frac{d}{dx} + \beta(x) $$
But here I'm stuck at Riccati which looks like
$$ -\beta'(x) + \beta(x)^2 = x^{2n} - \epsilon $$
I've tried a finite-degree polynomial, but it seems that there is no such solution for $\beta$ when $n>1$.
Finally, I've tried using Wolfram Mathematica to solve the initial equation, but it seems to be stuck after the $n=4$ (in cases of $n=2, 3$ it produces solutions in terms of Heun functions which I am unfamiliar with).