\[ -a x^3+x y'(x)+x y(x)^2-y(x)=0 \] ✓ Mathematica : cpu = 0.0203004 (sec), leaf count = 30
\[\left \{\left \{y(x)\to \sqrt {a} x \tanh \left (\frac {1}{2} \sqrt {a} \left (2 c_1+x^2\right )\right )\right \}\right \}\]
✓ Maple : cpu = 0.056 (sec), leaf count = 22
\[ \left \{ y \left ( x \right ) =\tanh \left ( {\frac {{x}^{2}+2\,{\it \_C1}}{2}\sqrt {a}} \right ) x\sqrt {a} \right \} \]
\begin {equation} xy^{\prime }+xy^{2}-y-ax^{3}=0\nonumber \end {equation} This is Riccati non-linear first order. But using the transformation \(y=xv\) it is transformed to easily solved ODE\[ y^{\prime }=v+xv^{\prime }\]
Therefore the ODE becomes
\begin {align*} x\left ( v+xv^{\prime }\right ) +x\left ( xv\right ) ^{2}-xv-ax^{3} & =0\\ xv+x^{2}v^{\prime }+x^{3}v^{2}-xv-ax^{3} & =0\\ x^{2}v^{\prime }+x^{3}v^{2}-ax^{3} & =0\\ v^{\prime }+xv^{2}-ax & =0\\ \frac {dv}{dx} & =x\left ( a-v^{2}\right ) \\ \frac {dv}{a-v^{2}} & =xdx \end {align*}
Integrating
\begin {align*} \frac {1}{\sqrt {a}}\tanh ^{-1}\left ( \frac {v}{\sqrt {a}}\right ) & =\frac {x^{2}}{2}+C\\ \tanh ^{-1}\left ( \frac {v}{\sqrt {a}}\right ) & =\sqrt {a}\left ( \frac {x^{2}}{2}+C\right ) \\ \frac {v}{\sqrt {a}} & =\tanh \left ( \sqrt {a}\left ( \frac {x^{2}}{2}+C\right ) \right ) \\ v & =\sqrt {a}\tanh \left ( \sqrt {a}\left ( \frac {x^{2}}{2}+C\right ) \right ) \end {align*}
Therefore
\begin {align*} y & =xv\\ & =x\sqrt {a}\tanh \left ( \sqrt {a}\left ( \frac {x^{2}}{2}+C\right ) \right ) \end {align*}
Verification
restart; ode:=x*diff(y(x),x)+x*y(x)^2-y(x)-a*x^3=0; my_solution:=x*sqrt(a)*tanh(sqrt(a)*(x^2/2+_C1)); odetest(y(x)=my_solution,ode); 0