Not sure if you have your heart set on exact solutions, but what about a perturbative solution? I can only speak in relatively general terms for now, but I hope this outline is clear enough to guide you on what to do.
Rewrite the first equation for which we know the solution as
$x y''+2 x y' + x^2 y = 0$
We write the desired solution to this equation as $y_0(x)$. You write a general form of the solution above; Mathematica, however, writes the form more compactly as
$y_0(x) = e^{-x} (a \, \mathrm{Ai}{[e^{i \pi/3} (1-x)]} + b \, \mathrm{Bi}{[e^{i \pi/3} (1-x)]}) $
$\mathrm{Ai}$ and $\mathrm{Bi}$ being the Airy functions and the constants $a$ and $b$ determined by initial conditions not specified.
We rewrite the 2nd equation as follows:
$x y''+2 x y' + x^2 y = (4 x+1) y'$
The first step in a perturbative approach is to introduce a "small" parameter $\epsilon$ that represents a deviation from the initial solution. That is, we write the the above equation as
$x y''+2 x y' + x^2 y = \epsilon (4 x+1) y'$
This is a regular perturbation because the perturbation is not being done at the highest derivative (i.e., no new solutions would be introduced by the perturbation). The solution $y(x)$ may then be written as
$y(x) = \sum_{n=0}^{\infty} \epsilon^n y_n(x)$
without really knowing whether that sum converges. The $y_n(x)$ satisfy a recurrence relation:
$x y_n''+2 x y_n' + x^2 y_n = \epsilon (4 x+1) y_{n-1}'$
where, for $n \ge 1$, $y_n(0) = 0$, $y_n'(0) = 0$.
Of course, you don't want to have to compute more than a few terms. You also need to determine whether the above series converges or diverges; if it diverges, then you need to figure out how many terms you can compute before the approximation no longer works. Finally, of course, to get an approximation to the solution, use $\epsilon = 1$.