So, I have a system of linear inhomogeneous PDE for a function $u(\rho, z)$: $$u_{,\rho} = \rho \phi_{,z}\\ u_{,z} = -\rho \phi_{,\rho}$$ with $\phi(\rho, z)$ being an axially symmetric solution of Laplace equation in cylindrical polar coordinates: $\phi_{,\rho\rho}+\frac{1}{\rho}\phi_{,\rho}+\phi_{,zz} = 0$.
My question is whether there is any procedure/trick to obtain a solution of this problem analytically.
For example, I've found out that it is possible to solve such an equation without any assumptions about $\phi(\rho, z)$: $$u_{,\rho}\phi_{,z}-u_{,z}\phi_{,\rho} = 0$$ Indeed, $u = Ф(\phi(\rho, z))$
I have tried reducing my problem to a single equation $u_{,\rho}\phi_{,\rho}+u_{,z}\phi_{,z} = 0$, but, unfortunately, do not see any sign of success in such a way. Method of characteristics won't work for general enough $\phi(\rho, z)$.
I will appreciate any kind of help or advice.