Solution (source code)

= Solution

On an interval avoiding zero put $v=y/x$. A real solution requires $|v|\leq1$, and the principal <square root> is $\sqrt{x^2-y^2}=|x|\sqrt{1-v^2}$. The <homogeneous first-order differential equation> reduces to
$$
v'=\frac{\sqrt{1-v^2}}{|x|}.
$$
Where $|v|<1$, <separation of variables> gives $\arcsin v=\operatorname{sgn}(x)\log|x|+C$. Thus the nonconstant pieces are
$$
\boxed{y=x\sin\bigl(\operatorname{sgn}(x)\log|x|+C\bigr),}
$$
restricted to intervals where the cosine of the phase is nonnegative. This sign condition is essential when substituting back into the original <square root>.

There are also the solutions $y=x$ and $y=-x$. They can join the nonconstant pieces with a continuous first <derivative>, because the transformed rate vanishes but is not <Lipschitz continuous> at $v=\pm1$. For example, on the positive axis every $A>0$ gives
$$
y(x)=\begin{cases}-x,&0<x\leq A,\\x\sin\bigl(\log(x/A)-\pi/2\bigr),&A<x<Ae^\pi,\\x,&x\geq Ae^\pi.\end{cases}
$$
The <sticking solutions of a bounded square-root differential equation> describe the same joining mechanism on either half-axis after using the increasing variable $t=\operatorname{sgn}(x)\log|x|$.