\[ \boxed { {x}^{2}{\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}y \left ( x \right ) -a \left ( \left ( y \left ( x \right ) \right ) ^{n}-y \left ( x \right ) \right ) =0} \]
Mathematica: cpu = 19.524479 (sec), leaf count = 27 \[ \text {DSolve}\left [x^2 y''(x)-a \left (y(x)^n-y(x)\right )=0,y(x),x\right ] \]
Maple: cpu = 0.702 (sec), leaf count = 67 \[ \left \{ y \left ( x \right ) ={\it ODESolStruc} \left ( {\it \_a},[ \left \{ {\frac {\rm d}{{\rm d}{\it \_a}}}{\it \_b} \left ( {\it \_a} \right ) = \left ( -{{\it \_a}}^{n}a+a{\it \_a} \right ) \left ( {\it \_b} \left ( {\it \_a} \right ) \right ) ^{3}- \left ( {\it \_b} \left ( { \it \_a} \right ) \right ) ^{2} \right \} , \left \{ {\it \_a}=y \left ( x \right ) ,{\it \_b} \left ( {\it \_a} \right ) ={\frac {1}{x{\frac {\rm d}{{\rm d}x}}y \left ( x \right ) }} \right \} , \left \{ x={{\rm e}^ {\int \!{\it \_b} \left ( {\it \_a} \right ) \,{\rm d}{\it \_a}+{\it \_C1}}},y \left ( x \right ) ={\it \_a} \right \} ] \right ) \right \} \]