\[ x^2+x y'(x)+y(x)^2=0 \] ✓ Mathematica : cpu = 0.0568126 (sec), leaf count = 32
DSolve[x^2 + y[x]^2 + x*Derivative[1][y][x] == 0,y[x],x]
\[\left \{\left \{y(x)\to \frac {x (-\operatorname {BesselY}(1,x)-c_1 \operatorname {BesselJ}(1,x))}{\operatorname {BesselY}(0,x)+c_1 \operatorname {BesselJ}(0,x)}\right \}\right \}\] ✓ Maple : cpu = 0.036 (sec), leaf count = 27
dsolve(x*diff(y(x),x)+y(x)^2+x^2 = 0,y(x))
\[y \left (x \right ) = -\frac {\left (c_{1} \operatorname {BesselY}\left (1, x\right )+\operatorname {BesselJ}\left (1, x\right )\right ) x}{c_{1} \operatorname {BesselY}\left (0, x\right )+\operatorname {BesselJ}\left (0, x\right )}\]
Hand solution
\[ 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