Boundary condition dependence on NDSolve inconsistent

Question

In attempting to solve a fourth order system I have encountered an issue with the way NDSolve uses boundary conditions. Consider the following three different attempts:

Nep = 1/2;
Nhp = 2;
Nen = 2;
Nhn = 1/2;
xmin = -2;
xmax = 2;

sol1 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[xmin] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol1], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

sol2 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[0] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol2], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

sol3 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[xmax] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol3], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

The difference between the three is just in the position at which a value of $\phi$ is specified. Given that there is no explicit dependence on $\phi$ in the system of equations, this shouldn't affect the solution up to the addition of a constant. However, only the second, where $\phi$ is specified at $x=0$ runs successfully. Why is this happening when all three problems are well-posed? Is there a fix?

Thanks in advance for any help.

Answer 1

UnitStep causes the problems. Substitute it by a very narrow ArcTan and it works fine. Like

sol1 = First@NDSolve[{2 \[Phi]''[x] == Ne[x] - Nh[x] - 3 (1/Pi ArcTan[x/10^-7]), -Nh[x] \[Phi]'[x] - Nh'[x] == 0, -Ne[x] \[Phi]'[x] + Ne'[x] == 0, \[Phi][xmin] == 0, \[Phi]'[xmin] == 0, \[Phi]'[xmax] == 0,  Ne[xmin] Nh[xmin] == 1}, {\[Phi], Ne, Nh}, 
                        {x, xmin, xmax}] 

May be you need higher WorkingPrecision.

Answer 2

I guess it's again the "Shooting" method that's not robust enough. The problem can be circumvented with the help of new-in-12 nonlinear FiniteElement method.

Since FiniteElement is chosen, the pre-processor of NDSolve can no longer handle the equation system (at least now), so let's help it a bit:

Nep = 1/2;
Nhp = 2;
Nen = 2;
Nhn = 1/2;
xmin = -2;
xmax = 2;

midsol=DSolve[{-Ne[x] ϕ'[x] + Ne'[x] == 0, -Nh[x] ϕ'[x] - Nh'[x] == 0}, {Nh, Ne},x][[1]]

csol = Solve[1 == Ne@xmin Nh@xmin /. midsol, C[1]][[1]]

neweq = 2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - UnitStep[-x] (Nep - Nhp) /. midsol /. 
   csol /. C[2] -> c
(*
2 ϕ''[x] == -c E^-ϕ[x] + E^ϕ[x]/c + (3 UnitStep[-x])/2 - (3 UnitStep[x])/2
 *)

psol = ParametricNDSolveValue[neweq, ϕ, {x, xmin, xmax}, c, 
  Method -> "FiniteElement"]

ϕ'[xmin] == 0, ϕ'[xmax] == 0 is omitted, because they're the default setting of FiniteElement in this case, you may check the document of NeumannValue for more info.

para1 = FindRoot[psol[c][xmin] == 0, {c, 1}]
(* {c -> 1.70937} *)

para2 = FindRoot[psol[c][0] == 0, {c, 1}]
(* {c -> 1.} *)

para3 = FindRoot[psol[c][xmax] == 0, {c, 1}]
(* {c -> 0.585012} *)

Plot[psol[c][x] /. {para1, para2, para3} // Evaluate, {x, xmin, xmax}]