# 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_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.

• Device Modeling? There are several monographs published by the Wien branch of Springer Verlag that describe comprehensively the theory and the finite difference schemes used to numerically solve the drift-diffusion equations coupled. with the nonlinear Poisson equation, for example "Semiconductor equations" authored by Markowich, Ringhofer and Schmeiser (1990) – Daniele Tampieri Feb 25 '20 at 18:47

## 1 Answer

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!

• As has already been pointed out, there are entire books about this type of problem. If your difficulties result from small diffusion coefficients or such, these resources certainly have hints how to deal with that. If you simply have a bug in your code, this is not the right place to get it debugged. – Michael Renardy Jul 24 '20 at 16:52