How can solve these system of differential equation? Assume that $U(\rho,z)$ and $V(\rho,z)$ functions are given. 
We have the system of equations for these functions in the form
$$
\frac{\partial^2 U}{\partial\rho^2}+\frac{1}{\rho}\frac{\partial U}{\partial\rho}+\frac{\partial^2 U}{\partial z^2} = 0 \ ,
$$
and
$$
\frac{\partial V}{\partial\rho} = \rho\left[\left(\frac{\partial U}{\partial\rho}\right)^2+\left(\frac{\partial U}{\partial z}\right)^2 \right]\ ,
\qquad \frac{\partial V}{\partial z} = 2\rho\frac{\partial U}{\partial\rho}\frac{\partial U}{\partial z}\ ,
$$
Can we use separation of the variable in the first equation? 
How about in the second and the third equations?
I will start the solution problem over here to show you exactly what I mean. The super-position principle, as you mentioned exists in linear systems. Hence, for the first part:

$ \frac{\partial^2 U}{\partial \rho^2}+\frac{ 1}{ \rho}\frac{\partial U}{\partial \rho}+\frac{\partial^2 U}{\partial z^2} = 0 $

with $U = R(\rho)Z(z)$

we get: 

$ R''+(1/\rho) R' = -\lambda^2 R$

$ R = a_1 Y_0(\lambda \rho) + a_2 J_0(\lambda \rho) $

Here, $J_0$ and $Y_0$ are Bessel functions of the first and second kind, respectively. 

$ Z'' = -\lambda^2 Z$

$ Z = b_1 \sin(\lambda z)+b_2 \cos(\lambda z)$

Next, as I tried to mention, $V$ depends only on $U$. There is no need of anything like the superposition principle:

$ \frac{\partial V}{\partial z} = 2\rho R'Z'(RZ)$

$ V = 2\rho R'R\big(\int Z'Z\partial z\big) +f(\rho) = \rho R'RZ^2 +f(\rho)$

$ \frac{\partial V}{\partial \rho} = \rho R''RZ^2 + R'RZ^2 + \rho (R')^2 Z^2 +f'(\rho)$

We also have: 
$ \frac{\partial V}{\partial \rho} = \rho ((R')^2Z^2+R^2(Z')^2)$

Then you can use the quantities $a_1,a_2,b_1,b_2,\lambda$ and $f(\rho)$ to make sense of the equality:

$\rho R''RZ^2 + R'RZ^2 + \rho (R')^2 Z^2 +f'(\rho) = \rho ((R')^2Z^2+R^2(Z')^2)$

$\rho (-R'/\rho -\lambda^2 R)RZ^2 + R'RZ^2 + \rho (R')^2 Z^2 +f'(\rho) = \rho ((R')^2Z^2+R^2(Z')^2)$

$ f'(\rho) = \rho R^2(\lambda^2 Z^2+(Z')^2)$

Probably you should check again, but I believe that $ (\lambda^2 Z^2+(Z')^2) $ is always constant (i.e. it does not depend on $z$). Therefore, you will get:

$ f'(\rho) = \rho R^2(\text{some constant in }b_1,b_2,\lambda)$

and you can solve for $V$
Snifkes, Thank you very much for your solution. Actually, I have solved this problem in the prolate coordinates, 
$$
\rho=\sqrt{(x^2-1)(1-y^2)} \ , z=xy
$$ 
and we get
$$
\left[\left(x^2-1\right)U_x\right]_x +
\left[\left(1-y^2\right)U_y\right]_y=0 \ ,
$$
$$
V_x = \frac{1-y^2}{x^2-y^2}\left[x(x^2-1)U_x^2-x(1-y^2)U_y^2
-2y(x^2-1)U_xU_y\right] \ ,
$$
$$
V_y = \frac{x^2-1}{x^2-y^2}\left[y(x^2-1)U_x^2-y(1-y^2)U_y^2
+2x(1-y^2)U_xU_y\right] \ ,
$$

Solution of the first equation is
$$
U(x,y)=\sum_{n=0}^{\infty}(a_n P_n(x)+b_n Q_n(x))(c_n P_n(x)+d_n Q_n(x))
$$
But my question was different: You are making separation of variable in the function $U(\rho,z)$, finally it depend on $\lambda$ which is separation constant and you need to sum over lambda like: $U= \sum_{\lambda} R_{\lambda}(\rho)Z_{\lambda}(z)$ - which is correct, but other two equations are not linear. As far as I know it is impossible to use summation over with separation constant for non linear equations. What do you think about it?