Solving non-linear ODE for divalent solution at a 1-D surface boudary I am trying to solve the following equation for a positively charged plane with charge density $\sigma$ at $z = 0$. 
$$ \phi''(z)=-\frac{e}{\epsilon \epsilon_0} \big(z_+n_{+}  e^{-\beta z_+ e\phi}-z_{-} n_{-} e^{\beta z_- e\phi} \big) $$
So far I have determined that the electric field and potential must go to zero at infinity, and that, assuming there is no electric field at $z < 0 $, I have the following boundary condition:
$$\phi'(0) = -\frac{\sigma}{\epsilon\epsilon_0}$$
I am not sure how to proceed with solving the equation given the two nonlinear terms.
 A: The Poisson-Boltzmann equation is principally valid when dealing with 1:1 electrolytes i.e. $z_- = -1$ and $z_+ = +1$.Then, as soon as divalent ions enter the system things become a bit more tricky and people (physicists at least) might jump at your throat if you use Poisson-Boltzmann in this case. I will therefore assume in the following that you do have a 1:1 electrolyte.
The usual thing is to use dimensionless variables for the field.
if you denote $\psi \equiv \beta e \phi$ then you get a rewriting of the PB equation as:
$\psi'' = \kappa^2 \sinh \psi \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:(1)$                    
with the boundary condition $\psi'(0) = -4\pi l_B \overline{\sigma}$
where $\overline{\sigma}$ is the surface charge density expressed in electric charge $e$ per unit length squared, $l_B \equiv \frac{e^2}{4\pi \epsilon \epsilon_0 k_B T}$ is the Bjerrum length expressed in SI units and $\kappa^2 = 8 \pi l_B n_b$ is the square of the inverse Debye length where $n_b = n_+ = n_-$ is your bulk salt concentration.
Before trying to solve the non linear case, it is important to notice that $\psi$ tells you the typical electrostatic potential energy of an ion with respect to $k_B T$. Basically, if $\psi$ is small everywhere, then your charged plate doesn't perturb very much the homogeneous distribution of ions or at least it does so in a linear way.
As a matter of fact, the above rewritten PB equation becomes $\psi'' \simeq \kappa^2 \psi$ in this case which is called the Debye-Huckel equation or the linearized PB equation. Its solution is then a trivial decaying exponential.
Let us go back to the full non linear problem for one plate.
The idea is to solve it by quadrature. To do so we first notice that if we multiply (1) by $\psi'$ then we get:
$\left(\frac{1}{2}(\psi')^2 - \kappa^2 \cosh \psi \right)' = 0 $
This implies that:
$\frac{1}{2}(\psi')^2 - \kappa^2 \cosh \psi = Cte = E $
I named this constant $E$ because if you are used to nonlinear mechanics then you call energy the function which, upon derivation with respect to the evolution variable, gives you back the equation of motion.
Let us assume there won't be any sign change in $\psi'$, we can then write (assuming the surface charge density is positive):
$\frac{d\psi}{dz} = -\sqrt{2E + 2\kappa^2 \cosh \psi}$
which leads to 
$dz = -\frac{d\psi}{\sqrt{2E + 2\kappa^2 \cosh \psi}}$ 
Upon integration we get:
$z = -\int_{\psi=\psi(0)}^{\psi(z)} \:\frac{d\psi}{\sqrt{2E + 2\kappa^2 \cosh \psi}}$ 
In the general case, $E$ can be anything and the solution has to be expressed as an elliptic integral (which is what is written above) with $E$ as a parameter that has be found by solving an implicit non linear equation derived from the boundary conditions. And therefore, even the simple problem of an electrolyte confined in between two plates is still an active field of research in mathematical physics and colloidal science.
In our case though, we implicitly chose the gauge potential so that both $\psi$ and $\psi'$ vanish at infinity. Because $E$ has a constant value throughout space, it means that it has to be $-\kappa^2$ everywhere since it is $-\kappa^2$ at infinity. We therefore get the simpler integral:
$z = -\int_{\psi=\psi_0}^{\psi(z)} \:\frac{d\psi}{\sqrt{2}\kappa \sqrt{-1+ \cosh \psi}} = \frac{1}{\kappa}\ln \coth \frac{\psi(z)}{4} - \frac{1}{\kappa}\ln \coth \frac{\psi(0)}{4}$.
Let us denote $\kappa z_0 = \ln \coth \frac{\psi(0)}{4}$ we therefore have:
$\kappa(z+z_0) = \ln \coth \frac{\psi(z)}{4} $
we now use the fact that $\coth \frac{x}{2} = \frac{1+e^{-x}}{1-e^{-x}}$ to rewrite the above equation as:
$e^{\kappa(z+z_0)} = \frac{1+e^{-\psi(z)/2}}{1-e^{-\psi(z)/2}}$
Doing the algebra step by step yields:
Step 1- $e^{\kappa(z+z_0)}(1-e^{-\psi(z)/2})=1+e^{-\psi(z)/2}$
Step 2- $-e^{-\psi(z)/2}(1+e^{\kappa(z+z_0)})=1-e^{\kappa(z+z_0)}$
Step 3- $e^{-\psi(z)/2} = \frac{e^{\kappa(z+z_0)}-1}{1+e^{\kappa(z+z_0)}}$
Step 4- $e^{-\psi(z)/2} =\frac{e^{-\kappa(z+z_0)}}{e^{-\kappa(z+z_0)}} \frac{e^{\kappa(z+z_0)}-1}{1+e^{\kappa(z+z_0)}} = \frac{1-e^{-\kappa(z+z_0)}}{1+e^{-\kappa(z+z_0)}}$
Final step- $\psi(z) = 2 \ln \frac{1+e^{-\kappa(z+z_0)}}{1-e^{-\kappa(z+z_0)}}$
Which is the general solution to the problem. The unique solution corresponding to the problem is found by looking for the value of $z_0$ that satisfies the boundary condition at the beginning.
