Solution (source code)

= Solution

Write $f(x,y)=F(u,v)$ in the proposed local coordinates. The <multivariable chain rule> gives
$$
f_x=2xF_u+2yF_v,\qquad f_y=-2yF_u+2xF_v.
$$
Consequently
$$
yf_x+xf_y=2(x^2+y^2)F_v=0.
$$
Away from the origin the <Jacobian determinant> is $u_xv_y-u_yv_x=4(x^2+y^2)>0$, so the <inverse function theorem> justifies these local coordinates. Therefore $F_v=0$ and the local <general solution> is
$$
\boxed{f(x,y)=\Phi(x^2-y^2).}
$$
Substitution verifies the result: $f_x=2x\Phi'$ and $f_y=-2y\Phi'$, whose weighted sum vanishes.

There is a global qualification. The coordinate map is two-to-one, and a <level set> of $x^2-y^2$ has separate <characteristic curve> branches. A single profile on the entire plane does not follow without a branch-identification assumption. For example, take $f=\exp[-1/(x^2-y^2)^2]$ on the wedge $x>0$, $x^2>y^2$, and zero elsewhere. This is smooth, satisfies the equation, and has different values at $(1,0)$ and $(-1,0)$ despite equal $u$. The <characteristic components of a hyperbolic first-order equation> explain this obstruction. At the origin the coordinate map is singular and the original equation reads $0=0$; the boxed answer is the intended local form on a regular branch, or a global form under an additional identification of branches.