Solving $ y\frac{\partial f}{\partial x} = x\frac{\partial f}{\partial y}$ rigorously I have the following PDE:
$$ y\frac{\partial f}{\partial x} = x\frac{\partial f}{\partial y} $$
For the sake of the problem, let us assume all functions in question are $C^{1}$-functions.
As far as as I understand, the most general solution is of the form $f(x, y) = g(x^{2} + y^{2})$ for any $C^{1}$-function $g:\mathbb{R}\rightarrow\mathbb{R}$.

My main question is, how can I prove this rigorously?


I have some work that would count as a non-rigorous derivation.
Start by separation of variables: $f(x, y) = X(x)Y(y)$. We get
$$ \frac{X'}{xX} = \frac{Y'}{yY} $$
The LHS is independent of $y$ and the RHS is independent of $x$, so both can be set to a constant $\lambda$.
This gives two ODEs
$$ X' - \lambda xX = 0 \qquad\text{ and }\qquad Y' - \lambda yY = 0. $$
These can be solved exactly with solutions
$$ X(x) = Ae^{-\lambda x^{2}/2} \qquad\text{ and }\qquad Y(y) = Be^{-\lambda y^{2}/2}. $$
Hence we can get a solution of the form
$$ f(x, y) = ABe^{-\lambda(x^{2} + y^{2})/2}, $$
or by "cleaning up the constants" we can write
$$ f(x, y) = C_{\lambda}e^{-\lambda(x^{2} + y^{2})}. $$
Now a general solution can be thought of as a "superposition" of uncountably many solutions of the above form, giving us a general solution
$$ f(x, y) = \int_{-\infty}^{\infty} C(\lambda)e^{-\lambda(x^{2} + y^{2})} \, d\lambda $$
for a general function $C:\mathbb{R}\rightarrow\mathbb{R}$.
This is the bilateral Laplace transform whose result is a general function $\widetilde{C}(x^{2}+y^{2})$.
I am wondering what would be a fully rigorous justification of all the steps here. I think the weakest part of my intuitive argument is the suggestion that our integral can return any function $\widetilde{C}(x^{2}+y^{2})$. I don't think this is even true. I'm sure there's also issues of convergence that I've handwaved away.
If there is a totally different way to prove this, I would be interested in seeing it. Any level of sophistication is ok. I just want to see what would be the most rigorous approach to this problem.
 A: This problem turns out to be solvable by a trick actually. There are some technicalities involving the differentiability class of $g$, but they are inessential to the main problem.
Suppose $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is a $C^{1}$-function such that $y\partial_{x}f = x\partial_{y}f$ (this is our PDE). For any given $r > 0$, define the path
$$ \gamma(t) = (r\cos t, r\sin t) $$
on $\mathbb{R}^{2}$. Then define $\sigma(t) = f\circ\gamma(t)$ as a function $\mathbb{R}\rightarrow\mathbb{R}$. This is a composition of $C^{1}$-functions, so it is itself a $C^{1}$-function. By the chain rule we have
$$ \sigma'(t) = \partial_{x}f(\gamma(t))\cdot -r\sin t + \partial_{y}f(\gamma(t))\cdot r\cos t = -y_{0}\partial_{x}f(x_{0}, y_{0}) + x_{0}\partial_{y}(x_{0}, y_{0}) = 0 $$
where the last equality comes from our PDE condition and $(x_{0}, y_{0}) = \gamma(t)$ for a given value of $t$. This is zero, so we infer that $\sigma(t)=\text{const}$. In particular, this can be easily (and rigorously) used to conclude that $f$ is constant on each set
$$ C(r) = \{ (x, y) : x^{2} + y^{2} = r^{2} \} $$
as long as $r > 0$. Since $C(0)$ is just the set consisting of the origin, $f$ is constant on that set as well. Thus, we may write $f(x, y) = g(x^{2} + y^{2})$ for some function $g:[0, \infty)\rightarrow\mathbb{R}$.
It only remains to show that $g$ is $C^{1}$. Unfortunately, this is not true for the entire domain (see the counterexample at the end). We can however show that $g$ is $C^{1}$ on $(0, \infty)$: We know $g(x) = f(\sqrt{x}, 0)$ for all $x > 0$, so $g$ restricted to $(0, \infty)$ is a composition of $C^{1}$-functions (on open subsets of Euclidean spaces) and so it is itself $C^{1}$ on $(0, \infty)$.
Counterexample. When I ask whether $g$ is $C^{1}$ on its entire domain $[0, \infty)$, I am asking whether it extends to a $C^{1}$-function on some $(-\varepsilon, \infty)$ with $\varepsilon>0$.
This is in general not possible, because of the $C^{1}$-function $f(x, y) = (x^{2}+y^{2})^{3/4}$ whose corresponding function here is $g(x) = x^{3/4}$. The first derivative is $g'(x) = 3/4\cdot x^{-1/4}\rightarrow\infty$ as $x\rightarrow 0+$, so $g$ does not extend in the desired way.
Note. I think $g$ would have been better behaved if we specified that $f(x, y) = g(\sqrt{x^{2}+y^{2}})$.
