# Boundary condition dependence on NDSolve inconsistent

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.

• In my run, your sol1 also successed by with a warning throwed out, sol3 failed... Nov 3, 2020 at 8:29
• I'm not really counting that as a success when it isn't giving the correct answer. Nov 3, 2020 at 8:43

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.

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}]