$$ (x^2+x)\frac{\mathrm d^2y}{\mathrm dx^2}+(2-x^2)\frac{\mathrm dy}{\mathrm dx}-(2+x)y=x(1+x)^2\tag*{}$$
This is the D.E. I'm trying to solve via the method of variation of parameter. As per this method, you're required to find the two linearly independent solutions, say, $y_1(x)$, $y_2(x)$ of the reduced D.E. $\displaystyle (x^2+x)\frac{\mathrm d^2y}{\mathrm dx^2}+(2-x^2)\frac{\mathrm dy}{\mathrm dx}-(2+x)y=0\tag*{}$
I found that $y_1(x)=1/x$ and $y_2(x)=e^x$ work. Their Wronskian: $\displaystyle W(1/x, e^{x}) = \begin{vmatrix}1/x & e^{x} \\ -1/x^2&e^{x}\end{vmatrix}=e^{x}\left(\frac{x+1}{x^2}\right)\tag*{}$
According to the formula derived in my notes, the general soln. is given by: $$y=Ay_1+By_2\\-y_1\int\frac{Xy_2}{W(y_1, y_2)}\mathrm dx+y_2\int\frac{Xy_1}{W(y_1, y_2)}\mathrm dx$$ where $X$ is the function on the RHS of the original D.E., in this case, it's $x(1+x)^2$.
$\displaystyle \begin{aligned} y & =\frac{A}{x} +Be^{x} -\frac{1}{x}\int \frac{x^{2} \cdotp x( 1+x)^{2} e^{x}}{e^{x}( 1+x)}\mathrm{d} x+e^{x}\int \frac{x^{2} \cdotp x( 1+x)^{2}}{e^{x}( 1+x) \cdotp x}\mathrm{d} x\\ & =\frac{A}{x} +Be^{x} -\frac{1}{x}\int x^{3}( 1+x)\mathrm{d} x+e^{x}\int e^{-x} x^{2}( 1+x)\mathrm{d} x\\ & =\frac{A}{x} +Be^{x} -\frac{1}{x}\left(\frac{x^{4}}{4} +\frac{x^{5}}{5}\right) -\left( x^{3} +4x^{2} +8x+8\right) \end{aligned} $
However, this soln. doesn't seem to match with Wolframalpha's:
I have double checked all my steps. I need help with this ambiguity. I am pretty confident that the formula isn't incorrect because I've used it to solve similar D.E.s and my answer checked with Wolframalpha every time.
The issue comes from your $X$ term. Before using Variation of Parameters method, you need to make the coefficient of the second derivative term to be $1$, namely,
$$\frac{\mathrm d^2y}{\mathrm dx^2}+\frac{(2-x^2)}{(x^2+x)}\frac{\mathrm dy}{\mathrm dx}-\frac{(2+x)}{(x^2+x)}y=\frac{x(1+x)^2}{(x^2+x)}=1+x$$
Hence, your $X$ term should be $$X=1+x$$