System of non-linear partial differential eqs from electrostatics

I have an electrostatics problem wich leads to the following system of differential equations:

$\displaystyle \frac{\partial E_z}{\partial z}=\frac{\rho}{\epsilon_0}$ **(1)**

$\displaystyle Z_i E_r \frac{\partial \rho}{\partial r}+(u_g+ Z_i E_z) \frac{\partial \rho}{\partial z} + \rho Z_i \frac{\partial E_z}{\partial z}=0$ ** (2)**

Substituting eq. (1) into eq. (2):

$\displaystyle Z_i E_r \frac{\partial \rho}{\partial r}+(u_g+ Z_i E_z) \frac{\partial \rho}{\partial z} + \frac{\rho^2 Z_i}{\epsilon_0}=0$ ** (3)**

Therefore I have a system of 2 equations (1 & 3) with 2 unknowns, the axial field $\displaystyle E_z$ and the charge density $\displaystyle \rho(z,r)$. The rest of the variables are known so they can be supposed as constants.

The Neumann boundary conditions in the axial axis, when r=0 are:

$\displaystyle \frac{\partial E_z}{\partial r}=0; \frac{\partial E_z}{\partial z} \neq 0$

$\displaystyle \frac{\partial \rho}{\partial r}=0; \frac{\partial \rho}{\partial z} \neq 0$

I'm not sure on how to solve it, I'm considering two options:

- derivate eq. (3) with respect to $\displaystyle z$ to substitute in eq. (1), but I donīt get rid of $\displaystyle E_z$ and the eq. (3) becomes more complicated.

- Solve by semi-implicit method, considering that $\displaystyle z=du_z/dt$, but since is an equation in partial derivatives I'm not sure on how to manage the term in $\displaystyle r$

I'm totally stuck on this, I'm asking for a direction of solving it, not for a solution, so any help would be grateful.

Thanks in advance.

P.S.: I forgot to say that the problem is solved in cylindrical coordinates (it occurs due to a discharge tip inside a cylindrical tube) and it has axial symmetry.