How to solve numerically a system of 3 interdependent non-linear ordinary differential equations?

Question

As per title, I need to solve this:

$$ \begin{cases} \frac{d^2V}{dx^2} = -\frac{q}{\epsilon}\left[p - n + \frac{N_0}{1+c_pp+c_nn}\right] \\\\ \frac{d}{dx}\left[\mu_nn\frac{dV}{dx} + D_n\frac{dn}{dx}\right] - \left[\frac{n-n_0}{\tau_n} - G_{op}\right] = 0 \\\\ \frac{d}{dx}\left[\mu_pp\frac{dV}{dx} - D_p\frac{dp}{dx}\right] + \left[\frac{p-p_0}{\tau_p} - G_{op}\right] = 0 \end{cases} $$

where V (electrostatic potential), n (total free electrons density) and p (total free holes density) are the unknown functions.

Being in 1D and in static conditions, all the functions depend on just the spatial variable "x" (the domain goes from "0" to "L").

Finally, the boundary conditions are:

$$ \begin{cases} V(x=0) = V_{left} (known) \\\\ V(x=L) = V_{right} (known) \\\\ \frac{dp}{dx}(x=0) = 0 \\\\ \frac{dn}{dx}(x=0) = 0 \\\\ \frac{dp}{dx}(x=L) = 0 \\\\ \frac{dn}{dx}(x=L) = 0 \end{cases} $$

It is basically a Poisson equation coupled with drift-diffusion + continuity equations for electrons and holes in a crystalline semiconductor.

Answer 1

I tried an iterative approach for which I set an initial guess for "n" and "p" (which are the functions I have at equilibrium), solve Poisson with boundary conditions slightly changed from equilibrium conditions, put the resulting potential profile in the 2 drift-diffusion + continuity equations and derive new guesses for holes and electrons densities to put again in the Poisson equation: iterate until it converges.

Then change again the potential boundary conditions and do it all again using as initial guess for "n" and "p" the results of the previous iterative procedure.

Finally, the potential boundary conditions are slowly brought to the desired value.

Practically, I solve the system at conditions slowly varying from equilibrium in an iterative way.

It converges but gives crazy (and wrong) results. I can provide the Matlab file!