# Drift-diffusion ODEs with two boundary conditions

I am trying to solve drift-diffusion equatons (Poisson's equation, continuity equations for electrons and holes and Kirchhoff's equation) for a reversely biased diode in a stationary state (no time dependences therefore). Recombination, impact ionization and drift velocity nonlinearity are taken into account. However, I need to implement four boundary conditions for ohmic contacts (two for each) instead of initial conditions. Commercial TCADs do it quite well, but I need Mathematica to solve it. My code is

q = 16*10^-20; t0 =
1 10^-6; A = 10^9; \[Epsilon] = 12; Subscript[\[Epsilon], 0] =
882*10^-16; W = 10^-2; R = 50; Vs = 10^7;
Esp = 23200;
Esn = 8000;
bn = 11*10^5;
bp = 22*10^5;
\[Alpha]ns = 74*10^4;
\[Alpha]ps = 725*10^3;
S = 10^-2; Er = 5/W;
Ur = 999;
\[Tau]n = 10*10^-6; \[Tau]p = 10*10^-6;
ni = 58*10^8;
T = 300;
k = 138*10^-25;
large = 100000;
\[Delta] = 10^-4;
SetSystemOptions["CheckMachineUnderflow" -> False];

Nd[x_] =
10^18*1/(1 + Exp[2*large*x]) - 10^18*1/(1 + Exp[2*large*(W - x)]);

Vn[x_] = Vs*F[x]/(Esn + F[x]); Vp[x_] = Vs*F[x]/(Esp + F[x]);
Difn[x_] = (k*T)/q*Vs*1/(Esn + F[x]);

Difp[x_] = (k*T)/q*Vs*1/(Esp + F[x]); \[Alpha]n[x_] = \[Alpha]ns*
Exp[-bn/(F[x])]; \[Alpha]p[x_] = \[Alpha]ps*Exp[-bp/(F[x])];

equ = {
0 == \[Alpha]n[x]*Vn[x]*n[x] + \[Alpha]p[x]*Vp[x]*p[x] +
D[Difp[x]*D[p[x], x], x] - (
p[x]*n[x] - ni^2)/(\[Tau]p*(n[x] + ni) + \[Tau]n*(p[x] + ni)) -
D[Vp[x]*p[x], x],
D[F[x], x] ==
q/(\[Epsilon]*Subscript[\[Epsilon], 0])*(Nd[x] + p[x] - n[x]), (
A*t0 - Ur)/(q*R*S) ==
Vp[x]*p[x] + Vn[x]*n[x] - Difp[x]*D[p[x], x] + Difn[x]*D[n[x], x],
n[0 - \[Delta]] == (Sqrt[(Nd[0 - \[Delta]])^2/4 + ni^2] +
Nd[0 - \[Delta]]/2),
p[0 - \[Delta]] ==
ni^2/(Sqrt[(Nd[0 - \[Delta]])^2/4 + ni^2] + Nd[0 - \[Delta]]/2),
n[W + \[Delta]] ==
ni^2/(Sqrt[(Nd[W + \[Delta]])^2/4 + ni^2] - Nd[W + \[Delta]]/2),
p[W + \[Delta]] == (Sqrt[(Nd[W + \[Delta]])^2/4 + ni^2] -
Nd[W + \[Delta]]/2)};

soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]}];


Which gives an error:

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

The same error occurs when I try shooting method.

When I do what Mathematica suggested

soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]},
Method -> {"EquationSimplification" -> "Residual"}];


the error becomes:

NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems.


Is there any way to solve my problem?

EDIT #1

When I use

soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]},
Method -> {"EquationSimplification" -> "Solve"}];


as xzczd suggested, I get a better result (at least something):

Power::infy: Infinite expression 1/0. encountered.

Infinity::indet: Indeterminate expression E^ComplexInfinity encountered.


Actually, the electric field F(x) must be positive everywhere, so I don't know what's wrong here.

• Don't follow that suggestion, just use Method -> {"EquationSimplification" -> "Solve"} or SolveDelayed->False to force NDSolve to find the explicit formula, NDSolve is capable of doing that in your case. Then you'll see some other warning, perhaps something is wrong with your model, perhaps you need to tackle the initial guess of Shooting method carefully, I haven't looked into it so I'm not sure, but at least you're in the correct direction now. Commented Apr 1, 2023 at 14:09
• Is it possible to normalize functions so that to exclude large parameters like 12500000000000000000 from equations? Commented Apr 1, 2023 at 15:53

Sorry, but again: Scale down your physical system. SI numbers on the atomic scale make no sense numerically:

In[23]:=
Nd[x_] =
10^18*1/(1 + Exp[2*large*x]) - 10^18*1/(1 + Exp[2*large*(W - x)])

Out[23]= -(1000000000000000000/(
1 + E^(200000 (1/100 - x)))) + 1000000000000000000/(1 + E^(200000 x))

In[25]:= -(1000000000000000000/(1 + E^(200000 (1/100 - x)))) +
1000000000000000000/(1 + E^(200000 x)) // N // Together

Out[25]= (1. (3.881180194284369*10^886 -
1.*10^18 2.71828^(400000. x)))/((1. + 2.71828^(
200000. x)) (3.881180194284369*10^868 + 1. 2.71828^(200000. x)))


Is a good idea to convert the input into the Atomic System by hand.

Best way of course is to use the Entity and Quantity system and shift the Quantity scaling eg to a diffusion constant of

10^-5 cm^2/s == 1  MyUnit.

• Thank you for your reply, but that did not help. Commented Apr 6, 2023 at 13:00
• Nothing else can help. its impossible for any machine to work with ((1 + E^(200000 (1/100 - x)))) + 1000000000000000000/(1 + E^(200000 x). Call this expr[x_]=%//N and then Plot[expr[x], {x, 0, 1}]. Nothing but an empty plot. All numerical procedures like Plot, NDSolve convert integers to reals and produce a lot of over- and underflow instances. Blank output is an easy way to handle such a potentially dangerous input. Commented Apr 6, 2023 at 17:29
• These exponents are just analytical approximation of Heavyside step function, where 200000 is just a sufficiently large number to make the transition between regions sharper. I used them instead of HeavisideTheta[1/100-x], which I find even less machine-friendly. And simple Plot function plots it quite well and smooth. Commented Apr 6, 2023 at 18:45
• By the way, to plot this function correctly you need to set the correct scale: Plot[Nd[x], {x, -0.01, 0.02}, PlotRange -> {{-0.001, 0.011}, {-2*10^18, 2*10^18}}] Commented Apr 6, 2023 at 19:04