Help solving differential equation in closed form - Damped Harmonic Oscillator I am attempting to solve analytically, a differential equation of the form
$$ -\alpha \frac{d^2y}{dx^2} + \beta y \frac{dy}{dx} + \gamma x^2 y = (\epsilon )y$$
Where, $\alpha $ , $ \beta$, $\gamma$ and $ \epsilon$ are constants.
The inclusion of the second term has thrown a spanner in the works. I attempted using Frobenius method but am unable to formulate a recursion relation. I have the solution of the equation with the second term excluded, and it yields a solution that depends on a Hermite polynomial.
I am hoping the solution will incorporate that too. Let me know if you have any resources you can
 A: I'm not sure if it's possible to obtain a "nice" solution. But if you're okay with approximation, you could use perturbation theory here, taking $\beta$ as a small parameter. Write $y$ as a combination of simpler functions:
$$
y(x) = y_0(x) + \beta y_1(x) + \beta^2 y_2(x) + \cdots = \sum_{k=0}^\infty \beta^ky_k(x).
$$
A priori, this sum is formal--there is no reason it should converge in any sense (unless of course the series terminates). Substitute into the original ODE:
$$
-\alpha\sum_{k=0}^\infty \beta^k\frac{d^2y_k}{dx^2} + \beta\left(\sum_{k=0}^\infty\beta^ky_k\right)\left(\sum_{k=0}^\infty\beta^k\frac{dy_k}{dx}\right)+ \gamma x^2\sum_{k=0}^\infty \beta^ky_k = \epsilon \sum_{k=0}^\infty \beta^ky_k.
$$
Equate coefficients according to powers of $\beta$. For example, looking at $\beta^0$ gives you the linear harmonic oscillator equation,
$$
-\alpha\frac{d^2y_0}{dx^2} + \gamma x^2y_0 = \epsilon y_0.
$$
You remarked that you already have the solution to this equation. Take that solution $y_0$ and examine the $\beta^1$ equation,
$$
-\alpha\beta\frac{d^2y_1}{dx^2} + \beta y_0\frac{dy_0}{dx} + \gamma x^2\beta y_1 = \epsilon\beta y_1 \implies -\alpha\frac{d^2y_1}{dx^2} + \gamma x^2y_1 = \epsilon y_1 + y_0\frac{dy_0}{dx}.
$$
We see that this equation is actually still linear in $y_1$ because we know $y_0$! In fact, it is just the same as the linear harmonic oscillator with the inhomogeneous term $y_0\frac{dy_0}{dx}$. In general, following this procedure gives that for the $\beta^k$ term, we have
$$
-\alpha\frac{d^2y_k}{dx^2} + \gamma x^2y_k = \epsilon y_k + f_k,
$$
where$f_k$ is a term depending only on $y_0, \cdots, y_{k-1}$ and their first derivatives (you could write a formula for $f_k$ using combinatorical methods). So now you can apply techniques of linear ODEs to solve the problem.
Again, there is no guarantee that this process terminates, or even that the sequence $y_k$ converges to the true solution $y$. But it may give you at least an idea of the asymptotic behavior of your true solution. Perhaps look at what the ladder operators do to the $f_k$.
