\[ x^2+x y'(x)+y(x)^2=0 \] ✓ Mathematica : cpu = 0.0163594 (sec), leaf count = 32
\[\left \{\left \{y(x)\to \frac {x \left (-c_1 J_1(x)-Y_1(x)\right )}{c_1 J_0(x)+Y_0(x)}\right \}\right \}\]
✓ Maple : cpu = 0.076 (sec), leaf count = 40
\[ \left \{ y \left ( x \right ) =-{\frac {{\it \_C1}\,x{{\sl Y}_{1}\left (x\right )}}{{\it \_C1}\,{{\sl Y}_{0}\left (x\right )}+{{\sl J}_{0}\left (x\right )}}}-{\frac {{{\sl J}_{1}\left (x\right )}x}{{\it \_C1}\,{{\sl Y}_{0}\left (x\right )}+{{\sl J}_{0}\left (x\right )}}} \right \} \]
\[ xy^{\prime }+y^{2}+x^{2}=0 \]
This is Riccati first order non-linear. Writing it in standard form and for \(x\neq 0\)\begin {align} y^{\prime } & =-x-\frac {1}{x}y^{2}\tag {1}\\ & =f_{0}+f_{1}y+f_{2}y^{2}\nonumber \end {align}
Where \(f_{0}=-x,f_{1}=0,f_{2}=-\frac {1}{x}\). Using standard substitution \(y=\frac {-u^{\prime }}{uf_{2}}\) changes the ODE to second order linear ODE
\begin {equation} y=\frac {xu^{\prime }}{u}\tag {2} \end {equation}
Hence
\[ y^{\prime }=\frac {u^{\prime }}{u}+x\frac {u^{\prime \prime }}{u}-\frac {x\left ( u^{\prime }\right ) ^{2}}{u^{2}}\]
Equating this to RHS of (1) gives
\begin {align*} \frac {u^{\prime }}{u}+x\frac {u^{\prime \prime }}{u}-\frac {x\left ( u^{\prime }\right ) ^{2}}{u^{2}} & =-x-\frac {1}{x}\left ( \frac {xu^{\prime }}{u}\right ) ^{2}\\ \frac {u^{\prime }}{u}+x\frac {u^{\prime \prime }}{u} & =-x\\ u^{\prime \prime }+\frac {1}{x}u^{\prime }+u & =0 \end {align*}
This is Lienard ODE. Since it is not constant coefficient ODE, the solution will be in Bessel functions, using Power series method. The solution is
\[ u=C_{1}\operatorname {BesselJ}\left ( 0,x\right ) +C_{2}\operatorname {BesselY}\left ( 0,x\right ) \]
But \(\frac {d}{dx}\operatorname {BesselJ}\left ( 0,x\right ) =-\operatorname {BesselJ}\left ( 1,x\right ) \) and \(\frac {d}{dx}\operatorname {BesselY}\left ( 0,x\right ) =-\operatorname {BesselY}\left ( 1,x\right ) \), hence
\[ u^{\prime }\left ( x\right ) =-C_{1}\operatorname {BesselJ}\left ( 1,x\right ) -C_{2}\operatorname {BesselY}\left ( 1,x\right ) \]
And from (2) the solution is
\begin {align*} y & =\frac {xu^{\prime }}{u}\\ & =x\frac {\left [ -C_{1}\operatorname {BesselJ}\left ( 1,x\right ) -C_{2}\operatorname {BesselY}\left ( 1,x\right ) \right ] }{C_{1}\operatorname {BesselJ}\left ( 0,x\right ) +C_{2}\operatorname {BesselY}\left ( 0,x\right ) }\\ & =-x\frac {C_{1}\operatorname {BesselJ}\left ( 1,x\right ) +C_{2}\operatorname {BesselY}\left ( 1,x\right ) }{C_{1}\operatorname {BesselJ}\left ( 0,x\right ) +C_{2}\operatorname {BesselY}\left ( 0,x\right ) } \end {align*}
Let \(C=\frac {C_{1}}{C_{2}}\) then
\[ y=-x\frac {C\operatorname {BesselJ}\left ( 1,x\right ) +\operatorname {BesselY}\left ( 1,x\right ) }{C\operatorname {BesselJ}\left ( 0,x\right ) +\operatorname {BesselY}\left ( 0,x\right ) }\]
Verification
restart; ode:=x*diff(y(x),x)+y(x)^2+x^2=0; my_sol:=-x*(_C1*BesselJ(1, x)+BesselY(1,x))/(_C1*BesselJ(0, x)+BesselY(0,x)); odetest(y(x)=my_sol,ode); 0