\[ x^3 (-y(x))+y'(x)+x y(x)^2-2 x=0 \] ✓ Mathematica : cpu = 0.375703 (sec), leaf count = 63
\[\left \{\left \{y(x)\to \frac {2 c_1 x^2+\sqrt {\pi } x^2 \text {erf}\left (\frac {x^2}{2}\right )+2 e^{-\frac {x^4}{4}}}{2 c_1+\sqrt {\pi } \text {erf}\left (\frac {x^2}{2}\right )}\right \}\right \}\]
✓ Maple : cpu = 0.291 (sec), leaf count = 51
\[ \left \{ y \left ( x \right ) ={\frac {1}{\sqrt {\pi }} \left ( \sqrt {\pi }{\it Erf} \left ( {\frac {{x}^{2}}{2}} \right ) {\it \_C1}\,{x}^{2}+{x}^{2}\sqrt {\pi }+2\,{{\rm e}^{-1/4\,{x}^{4}}}{\it \_C1} \right ) \left ( {\it Erf} \left ( {\frac {{x}^{2}}{2}} \right ) {\it \_C1}+1 \right ) ^{-1}} \right \} \]
\begin {align} y^{\prime }-yx^{3}+xy^{2}-2x & =0\nonumber \\ y^{\prime } & =2x+yx^{3}-xy^{2}\nonumber \\ & =P\left ( x\right ) +Q\left ( x\right ) y+R\left ( x\right ) y^{2}\tag {1} \end {align}
This is Riccati first order non-linear ODE with \(P\left ( x\right ) =2x,Q\left ( x\right ) =x^{3},R\left ( x\right ) =-x\). We can convert Riccati to Bernoulli which is easier to solve using the substitution \(u=x^{2}-y\) or \(y=x^{2}-u\)\begin {align*} u^{\prime } & =2x-y^{\prime }\\ & =2x-\left ( 2x+yx^{3}-xy^{2}\right ) \\ & =2x-\left ( 2x+\left ( x^{2}-u\right ) x^{3}-x\left ( x^{2}-u\right ) ^{2}\right ) \\ & =2x-\left ( 2x+\left ( x^{5}-ux^{3}\right ) -x\left ( x^{4}+u^{2}-2x^{2}u\right ) \right ) \\ u^{\prime } & =2x-\left ( 2x+\left ( x^{5}-ux^{3}\right ) -\left ( x^{5}+xu^{2}-2x^{3}u\right ) \right ) \\ & =2x-2x-\left ( x^{5}-ux^{3}\right ) +\left ( x^{5}+xu^{2}-2x^{3}u\right ) \\ & =-x^{5}+ux^{3}+x^{5}+xu^{2}-2x^{3}u\\ & =-ux^{3}+xu^{2} \end {align*}
This is of the form \(u^{\prime }=P\left ( x\right ) +Q\left ( x\right ) u+R\left ( x\right ) u^{2}\) and since \(P\left ( x\right ) =0\) then it is Bernoulli differential equation. (when \(P\left ( x\right ) \neq 0\) and \(R\left ( x\right ) \neq 0\) it is Riccati). To solve Bernoulli we always start by dividing by \(u^{2}\)\[ \frac {u^{\prime }}{u^{2}}=-\frac {1}{u}x^{3}+x \] Then we let \(\zeta =-\frac {1}{u}\), hence \(\zeta ^{\prime }=\frac {u^{\prime }}{u^{2}}\), therefore the above becomes\begin {align*} \zeta ^{\prime } & =x^{3}\zeta +x\\ \zeta ^{\prime }-x^{3}\zeta & =x \end {align*}
Integrating factor is \(e^{-\int x^{3}dx}=e^{-\frac {x^{4}}{4}}\), hence \[ d\left ( e^{-\frac {x^{4}}{4}}\zeta \right ) =xe^{-\frac {x^{4}}{4}}\] Integrating both sides gives\[ e^{-\frac {x^{4}}{4}}\zeta =\int xe^{-\frac {x^{4}}{4}}dx+C \] \(\int xe^{-\frac {x^{4}}{4}}dx=\frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) \), hence from above\begin {align*} e^{-\frac {x^{4}}{4}}\zeta & =\frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) +C\\ \zeta & =e^{\frac {x^{4}}{4}}\left ( \frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) +C\right ) \end {align*}
Since \(\zeta =-\frac {1}{u}\) then\[ u=-e^{-\frac {x^{4}}{4}}\left ( \frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) +C\right ) ^{-1}\]
And since \(y=x^{2}-u\) then\begin {align*} y & =x^{2}+e^{-\frac {x^{4}}{4}}\left ( \frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) +C\right ) ^{-1}\\ & =x^{2}+\frac {e^{-\frac {x^{4}}{4}}}{\frac {\sqrt {\pi }}{2}\operatorname {erf}\left ( \frac {x^{2}}{2}\right ) +C} \end {align*}
Verification
eq:=diff(y(x),x)+x*y(x)^2-x^3*y(x)-2*x = 0; sol:=x^2+ exp(-x^4/4)/(_C1+ sqrt(Pi)/2*erf(x^2/2)); odetest(y(x)=sol,eq); 0