System of non-linear partial differential eqs from electrostatics

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

**(1)**

** (2)**

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

** (3)**

Therefore I have a system of 2 equations (1 & 3) with 2 unknowns, the axial field and the charge density . 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:

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

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

- Solve by semi-implicit method, considering that , but since is an equation in partial derivatives I'm not sure on how to manage the term in

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.