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

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_nnfrac{dV}{dx} + D_nfrac{dn}{dx}right] – left[frac{n-n_0}{tau_n} – G_{op}right] = 0 \\
frac{d}{dx}left[mu_ppfrac{dV}{dx} – D_pfrac{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.