ndsolve solution goes to zero for large domain

I am trying to solve numerically the diffusion equation with absorbing boundary conditions on a finite domain. The initial condition is an approximation of a delta function centered in x0 with standard deviation sigma.

I find that changing the size of the domain slightly has a dramatic (non-physical) effect on the solution evaluated at a finite time. Why is this? Why don't I get a warning from Mathematica? This code reproduces the problem:

sigma = 1/32;
T = 10;
x0 = 1;
L = 24;
sol1 = NDSolve[{D[pn[x, t], {t, 1}] == 1/8 D[pn[x, t], {x, 2}],
pn[0, t] == 0, pn[L, t] == 0,
pn[x, 0] == 1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]},
pn[x, t], {x, 0, L}, {t, 0, T}];
L = 26;
sol2 = NDSolve[{D[pn[x, t], {t, 1}] == 1/8 D[pn[x, t], {x, 2}],
pn[0, t] == 0, pn[L, t] == 0,
pn[x, 0] == 1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]},
pn[x, t], {x, 0, L}, {t, 0, T}];
t = T;
GraphicsRow[{Plot[Evaluate[pn[x, t] /. sol1], {x, 0, 24},
PlotRange -> All],
Plot[Evaluate[pn[x, t] /. sol2], {x, 0, 26}, PlotRange -> All]}]

On the left I plot the solution with domain size L=24, on the right I plot the solution with L=26. As you can see, the change is dramatic (notice the y axis values) and Mathematica produces no warning.

How can I solve this issue and be able to solve the equation numerically for L larger than 24? • Tip; Setting the variable of integration t to a numeric value with t = T messes things up when you (or I) recompute the solutions. It's very inconvenient. You can avoid needing such a workaround by solving for pn instead of pn[x, t] in your NDSolve calls. Alternatively, you could use the code pn[x, t] /. sol /. t -> T. – Michael E2 May 16 '17 at 21:00

You have control the quality of the grid/mesh of the spatial variable. Either with "MaxStepSize":

NDSolve[{D[ppn[x, t], {t, 1}] == 1/8 D[ppn[x, t], {x, 2}],
ppn[0, t] == 0, ppn[L, t] == 0,
ppn[x, 0] ==
1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, ppn, {x, 0,
L}, {t, 0, T},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxStepSize" -> 1}}]

or with "MinPoints":

NDSolve[{D[ppn[x, t], {t, 1}] == 1/8 D[ppn[x, t], {x, 2}],
ppn[0, t] == 0, ppn[L, t] == 0,
ppn[x, 0] ==
1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, ppn, {x, 0,
L}, {t, 0, T},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> Ceiling@L}}]