Drift-diffusion ODEs with two boundary conditions

Question

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.

Answer 1

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.