# Neumann b.c are not satisfied in NDSolve PDE

I am solving a PDE (p[x,t]) with 1 spatial dimension {x,-10,10} and 1 time dimension {t,0,100}. There is an external function, U. The code reads

U[x_] := Sin[ 12*2 Pi/20 x];

sol = NDSolve[{
D[p[x, t], t] - D[p[x, t], x, x] + D[ U[x] p[x, t], x] - 1 +
p[x, t] == NeumannValue[0, x == -10] + NeumannValue[0, x == 10],
p[x, 0] == 1
}
, p, {t, 0, 100}, {x, -10, 10}][];


I want to impose dp/dx = 0 at x=-10 and x=10 at all times. Checking with:

state = NDSolveProcessEquations[{
D[p[x, t], t] - D[p[x, t], x, x] + D[ U[x] p[x, t], x] - 1 +
p[x, t] == NeumannValue[0, x == -10] + NeumannValue[0, x == 10],
p[x, 0] == 1
}
, p, {t, 0, 100}, {x, -10, 10}][];

state["FiniteElementData"][
"PDECoefficientData"]["DiffusionCoefficients"]
state["FiniteElementData"][
"PDECoefficientData"]["ConservativeConvectionCoefficients"]
state["FiniteElementData"][
(*
{{{{-1}}}}
{{{{0}}}}
{{{{0}}}}
*)


confirms that the Neumann value of 0 should only apply to the diffusive term, as the "ConservativeConvectionCoefficients" and "LoadDerivativeCoefficients" are 0. So the problem should be correctly imposing dp/dx = 0 at x=-10 and x=10. However

Table[D[(p /. sol)[x, t], x] /. {x -> -10}, {t, 0, 100}]
Table[D[(p /. sol)[x, t], x] /. {x -> 10}, {t, 0, 100}]


shows that the solution does not give a derivative of 0 for dp/dx at any of the timepoints for x=10,-10. Moreover the solution is just weirdly asymmetric, i.e. try:

Plot[(p /. sol)[x, 100], {x, -10, 10}]


How can I integrate this equation better and force the Neumann b.c to be satisfied at all the time points? Note: If I switch to:

U[x_] := Sin[2 Pi/20 x];


the dp/dx at the boundaries get closer to 0, like 10^-3, but this is not very good, and eventually, I want U to be more oscillatory than Sin[2 Pi/20 x].

You need to make the grid denser. Try e.g.

molfem[measure_ : Automatic] := {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> MaxCellMeasure -> measure}};

sol = NDSolve[{D[p[x, t], t] - D[p[x, t], x, x] + D[U[x] p[x, t], x] - 1 + p[x, t] == 0,
p[x, 0] == 1},
p, {t, 0, 100}, {x, -10, 10}, Method -> molfem[0.01]][];

Table[D[(p /. sol)[x, t], x] /. {x -> -10}, {t, 0, 100}]
(* {0., 2.78726*10^-7, 2.78725*10^-7, 2.78725*10^-7, 2.78725*10^-7, ... *)


The NeumannValue[0, …]` can be omitted, because they're the default setting.