求助：用有限元方法NeumannValue边界条件无法生效

二楼补充一下，代码中出现的一些函数和常数引用自之前的计算结果，本楼把相关代码一并贴出，需要先运行本楼的代码，再运行上述代码
<< NumericalCalculus`
Needs["NDSolve`FEM`"] ;
gauss[x_, x0_, \[Sigma]_] := Exp[-(x - x0)^2/(2 \[Sigma]^2)];
gauss[x_, x0_, \[Sigma]_] :=
Piecewise[{{Exp[-(x - x0)^2/(2 \[Sigma]^2)],
Abs[x - x0] <= 10 \[Sigma]}, {0, Abs[x - x0] > 10 \[Sigma]}}];
\[Epsilon] = Quantity[1, "ElectricConstant"];
e = Quantity[1, "ElementaryCharge"];
k = Quantity[1, "BoltzmannConstant"];
T = Quantity[300, "Kelvins"];
\[Epsilon]b = 7.1;
ni = Quantity[2 10^11, "Centimeters^-3"];
Wfi = Quantity[0.33, "eV"];
u = Wfi/(k T);
nb = ni E^u;
nnb = QuantityMagnitude[nb, "Centimeters^-3"]
Ldi = Sqrt[(\[Epsilon] \[Epsilon]b k T)/(
e^2 nb) ] (*先点一下下一行再点会本行，然后光标放在右侧方框处等待，才会跳出建议栏目*)
nLdi = QuantityMagnitude[Ldi, "Nanometers"];
nVs = -0.42;
Vs = Quantity[nVs, "Volts"](*单位是电压不是电子伏特*);
vs = Vs e/(k T)
F = Function[{x}, Sqrt[2 (Cosh[u + x]/Cosh[u] - x Tanh[u] - 1)]];
Pdensity = Quantity[2.09, "Millijoules/Centimeters^2"];
\[Lambda] = Quantity[618 3/2, "Nanometers"];
\[Alpha] = 1/\[Lambda];
R = 0.426;
\[Tau]L = Quantity[300, "Femtoseconds"];
n\[Tau]L = QuantityMagnitude[\[Tau]L, "Picoseconds"];
\[Mu] = Quantity[0.202, "Centimeters^2 /Volts Seconds"];
dif = k T \[Mu]/e;
\[Tau] = Ldi^2/dif;
n\[Tau] = QuantityMagnitude[\[Tau], "Picoseconds"]
hv = Quantity[1.55, "Electronvolts"];
耗尽区方程
vO = NDSolveValue[{v'[z] == F[v[z]], v[0] == vs}, v, {z, 0, 200},
WorkingPrecision -> 35, MaxSteps -> 10^6, AccuracyGoal -> 30]
n0[z_] := E^vO[z]
Plot[QuantityMagnitude[ni, "Centimeters^-3"] E^u n0[z/nLdi], {z, 0,
20 nLdi}, AxesLabel -> {"z(nm)", "n(cm^-3)"}, PlotRange -> All]
p0[z_] := E^(-2 u) E^(-vO[z])
Plot[QuantityMagnitude[ni, "Centimeters^-3"] E^u p0[z/nLdi], {z, 0,
20 nLdi}, AxesLabel -> {"z(nm)", "n(cm^-3)"}, PlotRange -> All]
eF0O = NDSolveValue[{D[eF0[z], z] == -n0[z] + p0[z] + 1 - E^(-2 u),
eF0[200] == 0}, eF0, {z, 0, 200}, WorkingPrecision -> 35,
MaxSteps -> 10^6, AccuracyGoal -> 30]
Plot[QuantityMagnitude[e nb Ldi/(\[Epsilon] \[Epsilon]b), "V/m"] eF0O[
z/nLdi], {z, 0, 20 nLdi}, AxesLabel -> {"z(nm)", "E(V/m)"},
PlotRange -> All]

……
\[CapitalGamma]2 = NeumannValue[0., x == 0 || x == 200]
你把 z 写成 x 了。顺便你这个NeumannValue的设置其实是多余的，因为有限元（FiniteElement）法默认使用为零的NeumannValue。
最大的问题其实是你那个 "PDEDiscretization" -> {Automatic, Automatic不能瞎用，不是什么选项的默认值都是Automatic！文档里对于Automatic在这种时候会起什么作用没解释，但是检查了一下，在PDEDiscretization里设置Automatica好像会导致后面的所有子选项被忽略。（使用s["Coordinates"] 检查，会发现只用了20个线元。）正确的做法是设置 "PDEDiscretization" -> {"MethodOfLines" ，还有，MaxCellMeasure别搞那么小，你设0.001等同于用了200/0.001==200000.个线元去解这个方程，MaxCellMeasure -> 0.1 差不多了。
PrecisionGoal也是多余的，不过这里没引发什么问题。
改完再画，会发现NeumannValue在初始时的误差略大，但是还是迅速趋向了零。你要实在不喜欢就把MaxCellMeasure再搞小点：
Plot[Derivative[0, 1][s][t, #], {t, -10 \[Tau]L/\[Tau], 5},
PlotRange -> All] & /@ {0, 200}