在初始时刻将一个波包放在rbs处考察波包随时间的演化
r1 = 0.05;
r2 = 5;
xx1 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r1;
xx2 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r2;
GLb = (1 + 45 rb^8 + 24 PEb rb^10 - 64 rb^9 TEb)/(
6 rb^4); rbs = 0.6804971882808477`;
rbl = 2.780560197330779`; PEb = 0.2379247924446251`; TEb = 0.5`;
pdeSoln =
NDSolveValue[{D[\[Rho][t, rb], t] ==
500 TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb],
rb], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbs)^2/0.01)), xx1, xx2}, \[Rho][t,
rb], {t, 0, 20000}, {rb, r1, r2}]
Plot3D[pdeSoln, {t, 0, 20000}, {rb, r1, r2}, PlotRange -> {0, 6},
AspectRatio -> 10/13, PlotPoints -> 150, MaxRecursion -> 2,
ColorFunction -> "Rainbow"]
给出了一个有问题的结果

GLb图如下所示，rbs和rbl是对应两个GLb极小值

如果把波包放在rbl处，即\[Rho][0, rb] == 1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbl)^2/0.01))
就会得到好结果，请问问题出在哪

1. 你为什么觉得这个结果是有问题的？你所期待的结果是啥？
2. 你确定方程本身正确吗？可否给一下出处？
3. 默认情况下 eerr 警告的误差显然大了些，而尝试加密空间网格后时间步进的速度开始变得很小。（在我写这答案这段时间里还没算完。）你确定时间区间的设置合理吗？
4. 你是不是有一层楼被吞/仅自己可见了？是的话最好申诉下。

方程是Fokker planck方程，ρ是演化波包，波包有一个振幅，我选的振幅是1/（0.1√π），即小于6，但结果演化到后面发现波包振幅随着时间的增长都越来越大

1.我所期待的结果是3楼的结果两个峰交换位置，演化时间可能不一样。
2.方程是Fokker Planck方程，两个边界方程，一个初始波包。GLb是势函数。
3.时间区间我不确定。
4.3楼结果已显示

5. 你的两个边界是否是对无穷远的近似？

3楼代码是
r1 = 0.05;
r2 = 5;
xx1 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r1;
xx2 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r2;
pdeSoln =
NDSolveValue[{D[\[Rho][t, rb], t] ==
500*TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb],
rb], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbl)^2/0.01)), xx1, xx2}, \[Rho][t,
rb], {t, 0, 20000}, {rb, r1, r2}]
Plot3D[pdeSoln, {t, 0, 20000}, {rb, r1, r2}, PlotRange -> {0, 6},
AspectRatio -> 10/13, PlotPoints -> 150, MaxRecursion -> 2,
ColorFunction -> "Rainbow"]

……折腾了半天发现还是网格不够密。添加选项
Method -> mol[2000, 4]
即可，其中mol的定义是
mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
2 GHz的电脑，耗时约 60 秒；边界其实可以用两个为 0 的第一类边界代替，求解速度略快，50 秒。（不过我只测了一次。）

回复 xzcyr :两个波包演化到最后都应该达到各自的稳定值，但无论是2000格还是3000格都有一个波包达不到

网格2000和3000的图像，最后都没有达到平稳值

……坑了个爹的，NDSolve好像又出现了性能退化。同样的代码（mol[3000, 4]），版本9画出来是这样的：
Plot[pdeSoln /. {{rb -> rbs}, {rb -> rbl}} // Evaluate, {t, 0, tend}]
pdeSoln[[0]]["Coordinates"][[1]] // ListPlot
版本12.2是这样的：
12.2明显不对劲。

嗬，这个问题居然会用到这个我是真没想到。
仔细试验后发现，本问题若想高效求解，必须使用特殊的方式进行空间离散，此问题在Stackexchange帖子《How to solve the tsunami model and animate the shallow water wave?》（编号107988）里有所论及，简单地说就是，对于某些基于守恒率的偏微分方程，如果空间离散不将守恒率以某种方式纳入考量，数值求解就会出问题。（求解效率低下，甚至无法得出正确解。）以往见到的例子，守恒率内部的项貌似都具有某种非线性，可LZ的这个方程只是变系数，居然也有类似问题，还真是让我开了眼了。
总之上代码。离散将使用 pdetoode 函数（mathematica.stackexchange.com/a/127997/1871）进行辅助。没改的参数我就不重复贴了：
{eq, ic, bc} = {D[ρ[t, rb], t] == 500 TEb D[mid[t, rb], rb], ρ[0, rb] ==
1/(0.1 Sqrt[π]) E^(-((rb - rbs)^2/0.01)), ρ[t, rb] ==
0 /. {{rb -> r1}, {rb -> r2}}};
eqmid = mid[t, rb] == E^(-GLb/TEb) D[E^(GLb/TEb) ρ[t, rb], rb];
(*eqmid2=mid2[t,rb]==E^(GLb/TEb) ρ[t,rb];*)
tend = 2 10^4;
domain = {r1, r2};
points = 3000; difforder = 2;
grid = Array[# &, points, domain];
ptoofunc = pdetoode[{ρ, mid}[t, rb], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = Block[{mid}, Evaluate@ptoofunc[eqmid[[1]]] = ptoofunc[eqmid[[2]]];
ptoofunc[eq] // del];
odeic = ptoofunc[ic];
odebc = With[{sf = 1}, Map[sf
# + D[#
, t] &, bc, {2}] // ptoofunc];
sollst = NDSolveValue[{ode, odeic, odebc}, ρ /@ grid, {t, 0, tend}, jianshi[t](*,
SolveDelayed\[Rule]True*)]; // AbsoluteTiming
sol = rebuild[sollst, grid]
计算只需约15 s。结果：
Plot[sol[t, rb] /. {{rb -> rbs}, {rb -> rbl}} // Evaluate, {t, 0, tend}]
Plot3D[sol[t, rb], {t, 0, tend}, {rb, r1, r2}, PlotRange -> All,
PlotPoints -> 50]
可算像话了嗯。

……我傻了。仔细想想，有限元法的标准形式就是基于守恒率的啊：
bc = \[Rho][t, rb] == 0 /. {{rb -> r1}, {rb -> r2}};
pdeSoln = NDSolveValue[{D[\[Rho][t, rb], t] ==
500 TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb], rb], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbs)^2/0.01)), bc}, \[Rho][t, rb], {t, 0,
tend}, {rb, r1, r2},
Method -> {MethodOfLines,
SpatialDiscretization -> {FiniteElement,
MeshOptions -> MaxCellMeasure -> (r2 - r1)/3000}}]; // AbsoluteTiming
耗时只要1秒不到。结果和13楼是一样的。

@xzcyr 之所以不令边界处波包等于零，是因为GLb函数一旦变化的话，边界处相对来讲可能并不趋于0。另外此处波包指的是概率密度，直接令概率密度等于0，意味着吸收，我的边界条件是概率流密度为0，即反射。