For this second-order linear ODE (eigenvalue `e`

) $ $ y”(x) – 2g\, y'(x) + [-e + g^2 – (\frac{b^2x^2}{2}+a)^2 + b^2x]\, y(x)=0$ $ I know the correct analytical general solution form with two undetermined constants `c1,c2`

(verified in the code below where `y1,y2`

is given) $ $ y=c1\,y1(x)+c2\,y2(x)$ $ I want to impose the boundary condition $ y(\pm5)=0$ in order to get the eigenvalue $ \lambda$ . (Actually $ y(\pm\infty)=0$ , but not much numerical difference, I suppose.) So I put `c1=1`

and tried `FindRoot`

`e`

and `c2`

together, but it didn’t work well. So I was wondering if I need to tune any option or some other methods?

I also add a direct numerical solution as a reference.

`F := (D[#, {x, 2}] - 2 g D[#, x] + (-e + g^2 - (b^2 x^2/2 + a)^2 + b^2 x) #) &; y1[x_] := E^(-(1/6) b^2 x^3 - x a + x g) HeunT[e, 0, -2 a, 0, -b^2, x]; y2[x_] := E^(-(1/6) b^2 x^3 - x a + x g) HeunT[e, 0, -2 a, 0, -b^2, x]; y[x_] := c1 y1[x] + c2 y2[x]; F[y[x]] == 0 // FullSimplify g = 0.2; a = 0.3 - I 0.3; b = 1; FindRoot[0 == y[x] /. {{c1 -> 1, x -> 5}, {c1 -> 1, x -> -5}}, {{e, -0.3}, {c2, 1}}, WorkingPrecision -> 55] FF = F /. e -> 0; {vals, funs} = NDEigensystem[FF[yy[x]], yy[x], {x, -5, 5}, 5]; vals Plot[Evaluate@Abs@funs, {x, -5, 5}, PlotRange -> All, PlotLegends -> Automatic] `