\[ xyy^{\prime }=y^{2}+x\sqrt {4x^{2}+y^{2}}\] Solving for \(y\) gives\begin {align*} {\normalsize y} & {\normalsize =}\operatorname {RootOf}\left ( \_z^{4}-4+\left ( p^{2}-1\right ) \_z^{2}-2\_z^{3}p\right ) {\normalsize x}\\ y & =xf+g \end {align*}
Where \(f=\operatorname {RootOf}\left ( \_z^{4}-4+\left ( p^{2}-1\right ) \_z^{2}-2\_z^{3}p\right ) \) and \(g=0\). Taking derivative of the above w.r.t. \(x\) gives\begin {align*} p & =\left ( f+xf^{\prime }\frac {dp}{dx}\right ) +\left ( g^{\prime }\frac {dp}{dx}\right ) \\ p & =f+xf^{\prime }\frac {dp}{dx}\\ p-f & =xf^{\prime }\frac {dp}{dx} \end {align*}
Using values for \(f\ \)the above simplifies to \begin {equation} p-\operatorname {RootOf}\left ( \_z^{4}-4+\left ( p^{2}-1\right ) \_z^{2}-2\_z^{3}p\right ) =\left ( x\frac {d}{dp}\operatorname {RootOf}\left ( \_z^{4}-4+\left ( p^{2}-1\right ) \_z^{2}-2\_z^{3}p\right ) \right ) \frac {dp}{dx} \tag {2A} \end {equation} The singular solution is found by setting \(\frac {dp}{dx}=0\) which results in \(p=\operatorname {RootOf}\left ( \_z^{4}-4+\left ( p^{2}-1\right ) \_z^{2}-2\_z^{3}p\right ) \). Substituting this in (1) does not generate any real solutions (only 2 complex ones) hence will not be used.
The general solution is found by finding \(p\) from (2A). Since (2A) is not linear in \(p\), then inversion is needed. Writing (2A) as\begin {align*} \frac {dx}{dp} & =\frac {xf}{p-f}\\ \frac {1}{x}dx & =\frac {f}{p-f}dp \end {align*}
Due to complexity of result, one now needs to obtain explicit result for \(\operatorname {RootOf}\) which makes the computation very complicated. So this is not practical to solve by hand. Will stop here. It is much easier to solve this ode as a homogeneous ode instead which gives the solution as\[ -\frac {\sqrt {4x^{2}+y^{2}}}{x}+\ln \left ( x\right ) =c_{1}\]