I believe that a catenoid supports a parametrization $\sigma : U \subset \mathbb{R}^2 \rightarrow \mathbb{R}^3$ that forms a conjugate system (i.e.\, $\sigma_{uv} \in\mathrm{span}(\sigma_u, \sigma_v)$) with the additional property that the coordinate curves are geodesics. I tried the following computational approach in building such a prameterization. I chose the following usual parameterization of a catenoid:

$$ \sigma = \left(\begin{array}{c} \cosh(u) \cos(v),\\
                                    \cosh(u) \sin(v),\\
                                    u\end{array}\right) $$  

I assumed that there is a coordinate transformation $u = \phi(x,y)$ and $ v = \psi(x,y)$ such that they would give me the result I need.
Upon computing the following conditions that ensure the geodesic-conjugate properties

 - conjugate property $\det(\sigma_x,\sigma_y,\sigma_{xy}) = 0$
 - geodesic x-coordinates $\det(N,\sigma_x,\sigma_{xx}) = 0$
 - geodesic y-coordinates $\det(N,\sigma_y,\sigma_{yy}) = 0$

where $N$ is the Gauss map, I arrived to the following system of PDEs:

 - $\ln \left(\psi_{x}\right)-\ln \left(\phi_{x}\right) = 
\ln \left(\psi_{y}\right)-\ln \left(\phi_{y}\right)
$
 - $\frac{{\partial}}{{\partial}y}\ln \left(\psi_{y}\right)-\frac{{\partial}}{{\partial}y}\ln \left(\phi_{y}\right)
 = 
\frac{\left(\phi_{y}^{2}+\psi_{y}^{2}\right) \tanh \left(\phi \right)}{\phi_{y}}$

 - $\frac{{\partial}}{{\partial}x}\ln \left(\psi_{x}\right)-\frac{{\partial}}{{\partial}x}\ln \left(\phi_{x}\right)=\frac{\left(\phi_{x}^{2}+\psi_{x}^{2}\right) \tanh\left(\phi \right)}{\phi_{x}}$

Can someone clarify how I can continue from here? Unfortunately I don't know much about PDEs but I still believe that there should be an answer to this.