问题如图
下面是我计算的代码：
DiffusionEquation={D[\[Rho][r,t],t]-d*1/r*D[D[\[Rho][r,t],r]*r,r]==0,\[Rho][r,0]==c,Piecewise[{{D[\[Rho][r,t],r]==c,r<R0},{D[\[Rho][r,t],r]==0,r>R0}}]};
DSolve[DiffusionEquation,\[Rho][r,t],{r,0,+\[Infinity]},{t,0,+\[Infinity]}]
这样的边界条件应该怎么处理？求助各位！

首先，LZ对分段的初始条件的定义方法显然是错误的，正确的做法恰恰是，就按你图上的来：
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
然后，目前DSolve在处理无限域问题时的语法并不是在自变量定义域的位置写Infinity，而是，就直接空着。这部分内容参看DSolve的自带帮助的“范围”部分，这里不展开讲。因为——上面这两点问题相较于这个问题的真正麻烦之处根本就是微不足道。
这个问题的真正麻烦之处在于：
1. DSolve在经历了10.3版的增强之后，似乎依旧不能直接求解这个问题。
2. 极坐标无限域下的Laplacian其实可以使用汉克尔变换消掉微分，但是遗憾的是目前Mathematica的HankelTransform性能还很糟，当然了，这个还不算太难变通。
3. 最终我们所得的变换后的解，不知道是太复杂了怎么的，总之，倒不回来。当然了，我们依旧得到了一个含了积分号的解析解。
好的，我们开始求解。首先利用拉普拉斯变换消掉关于t的导数项：
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
teqn = LaplaceTransform[{DiffusionEquation[[1]]}, t, s] /.
Rule @@ DiffusionEquation[[2]] /.
HoldPattern@LaplaceTransform[a_, __] :> a /. \[Rho] -> (rho[#] &)
(*
这里做了一些字母替换，主要是为了方便阅读。大家记住rho[r]代表\[Rho][r, t]对t的拉普拉斯变换就行了。然后呢，假如HankelTransform具备和LaplaceTransform相当的性能的话，我们只要使用类似上面的步骤，就可以消掉关于r的导数了，但是，遗憾的是：
*)
-d HankelTransform[rho'[r]/r + rho''[r], r, w]
HankelTransform[-(rho'[r]/r + rho''[r]), r, w]
(*嗯你没看错，只是加了个负号，HankelTransform就跪了。更囧的是：*)
Assuming[{R0 > 0}, HankelTransform[Piecewise[{{c, r < R0}}, 0] , r, w]]
(*返回未计算表达式*)
(*然而这个算式和自带帮助里的例子其实是差不多的：*)
(*并且Integrate其实也是会算这个积分的：*)
term = Assuming[{R0 > 0},
Integrate[r Piecewise[{{c, r < R0}}, 0] BesselJ[0, r w], {r, 0, Infinity}]]
(* (c R0 BesselJ[1, R0 w])/w *)
(*
总之呢，经历以上步骤我们总算是手动拼出了Hankel变换后的方程。我们把它解出来：
*)
ttsol = p /. First@Solve[-term + s p + d w^2 p == 0, p]
(* (c R0 BesselJ[1, R0 w])/(w (s + d w^2)) *)
(*然后拉普拉斯逆变换：*)
tsol = InverseLaplaceTransform[ttsol, s, t]
(*
(c E^(-d t w^2) R0 BesselJ[1, R0 w])/w
*)
(*
现在，我们来到了最后的难关，汉克尔逆变换。很遗憾，这个符号解我们算不出来：
*)
InverseHankelTransform[tsol, w, r]
Block[{d = 1, c = 1, R0 = 1, t = 1}, InverseHankelTransform[tsol, w, r]]
Block[{d = 1, c = 1, R0 = 1, t = 1}, Integrate[w tsol BesselJ[0, w r], {w, 0, Infinity}]]
(*以上方法全部失效*)
(*不过我们总算是得到了一个含有积分号的解析解：*)
(*这个解析解当然也是可以用的。我们和数值解对比一下：*)
asol[t_, r_, c_, d_, R0_?NumericQ] :=
NIntegrate[w (c E^(-d t w^2) R0 BesselJ[1, R0 w])/w BesselJ[0, w r], {w, 0, Infinity}]
aplot = Plot[asol[1, r, 1, 1, 1], {r, 0, 3}, PlotRange -> All]; // AbsoluteTiming
(* {6.04077, Null} *)
nsol = Block[{d = 1, c = 1, R0 = 1},
NDSolveValue[DiffusionEquation, \[Rho], {t, 0, 2}, {r, 0, 10},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.01}}]]
(*注意，当有限元方法被用于NDSolve的空间离散时，如果不明确指定边界，NDSolve会自动使用等于0的NeumannValue做为边界——我们在此需要的就是这种边界。*)
(* 对比一下： *)
aplot~Show~Plot[nsol[r, 1], {r, -3, 3}, PlotStyle -> {Red, Dashed}]

DiffusionEquation = {D[\[Rho][r, t],
t] - (d D[D[\[Rho][r, t], r] r, r])/r == 0, \[Rho][r, 0] ==
Piecewise[{{c, r < R0}}]};
DSolve[DiffusionEquation, \[Rho], {r, t}]
14版本可以解析解了，如果无穷上限积分能用的话

近期重看了一点格林函数相关的东西（以前上学的时候就没整明白），发现这个问题用格林函数法来求解也不赖。虽然说如6楼所述，现在DSolve已经能直接解这个问题了，这里还是提一下。
首先计算与这个问题对应的基本解。基本解对应的初值问题是：
(*是的，我有意用了完整的二维坐标，因为实践表明，在直角系下处理这个问题更容易——汉克尔变换的麻烦之处3楼已有解释。*)
sys = With[{ρ = ρ[x, y, t]}, {D[ρ, t] -
d Laplacian[ρ, {x, y}] == 0, ρ == DiracDelta[x, y] /.
t -> 0}]
(*
GreenFunction的性能目前还比较拉，截止14.3版，它只可以直接计算非齐次项为 DiracDelta[x-y] DiracDelta[t- τ] 的二维热传导方程的基本解——结合人工分析把所需的基本解弄出来或许也行，但是为了尽可能减少求解中的人工分析，这里将使用傅立叶变换求解这个问题——FourierTransform当前的实现其实也挺拉的，但是给它打补丁比较容易，我们为它写个壳：
*)
Clear@ft;
ft[(h : List | Plus | Equal)[a__], t_, w_] := ft[#, t, w] & /@ h[a]
ft[a_ b_, t_, w_] /; FreeQ[b, Alternatives @@ t] := b ft[a, t, w]
ft[HoldPattern@
Derivative[a___, b_ /; b > 0, c___][u_][L___, t_, R___], {tL___,
t_, tR___}, {wL___, w_, wR___}] /;
Length@{a} == Length@{L} && Length@{tL} == Length@{wL} := (-I w)^
b ft[Derivative[a, 0, c][u][L, t, R], t, w]
ft[a_, t_, w_] := FourierTransform[a, t, w]
(* 然后就可以把这个偏微分方程化成常微分方程了： *)
tsys = ft[sys, {x, y}, {kx, ky}] /.
HoldPattern@
FourierTransform[a_, __] :> (a /. ρ -> (Ρ[#3] &))
(*
{-d (-kx^2 Ρ[t] - ky^2 Ρ[t]) + Ρ'[
t] == 0, Ρ[0] == 1/(2 π)}
*)
(*常微分方程很好办，直接DSolve：*)
tsol = DSolveValue[tsys, Ρ[t], t]
(* E^(-d kx^2 t - d ky^2 t)/(2 π) *)
(* 等一下，既然现在DSolve本来就可以直接解这个偏微分方程，我们这里是不是避免使用DSolve比较好？行，我们还是用拉普拉斯变换吧： *)
tsol = InverseLaplaceTransform[
SolveValues[LaplaceTransform[tsys[[1]], t, s] /. Rule @@ tsys[[2]],
LaplaceTransform[Ρ[t], t, s]][[1]], s, t]
(* E^(-d (kx^2 + ky^2) t)/(2 π) *)
(* 最后再用傅立叶逆变换变回去： *)
solfunc[r_] =
Simplify[
InverseFourierTransform[tsol, {kx, ky}, {x, y}], {t > 0,
d > 0}] /. x^2 + y^2 -> r^2
(* E^(-(r^2/(4 d t)))/(4 d π t) *)
(* 最后一步再对初始条件积分： *)
dis =
Sqrt[# . #] &[
Subtract @@ FromPolarCoordinates /@ {{r, 0}, {ri, θi}}] //
Simplify
(* Sqrt[r^2 + ri^2 - 2 r ri Cos[θi]]
sol = Integrate[ri solfunc[dis] c, {ri, 0, R0}, {θi, -Pi, Pi}] *)
(* Integrate[(c*ri*BesselI[0, (r*ri)/(2*d*t)])/
(E^((r^2 + ri^2)/(4*d*t))*(2*d*t)), {ri, 0, R0}] *)
(* 依旧是个含了积分的解，但是它在有限域上，（大概）比上面那个无限域的好算一点。 *)
(*最后照例验算一下：*)
nsol = Block[{d = 1, c = 1, R0 = 1},
NDSolveValue[DiffusionEquation, ρ, {t, 0, 2}, {r, 0, 10},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {FiniteElement,
MeshOptions -> MaxCellMeasure -> 0.01}}]];
Block[{d = 1, c = 1, R0 = 1,
Integrate =
NIntegrate[##, Method -> {Automatic, "SymbolicProcessing" -> 0}] &,
t = 1},
Plot[{sol, nsol[r, t]}, {r, 0, 3}, PlotStyle -> {Automatic, Dashed}]]

还是说一下怎么从GreenFunction的输出里把所需的基本解搞出来。DiracDelta[x-xi]DiracDelta[y-yi] DiracDelta[t- ti]这个非齐次项其实等效于把初始条件DiracDelta[x-xi]DiracDelta[y-yi]给嵌到了ti时刻，所以用下面的代码也可以得到基本解：
Simplify[
GreenFunction[D[u[x, y, t], t] - d Laplacian[u[x, y, t], {x, y}],
u[x, y, t], {x, y} ∈ FullRegion[2](*{{x,-Infinity,
Infinity},{y,-Infinity,Infinity}}居然不行*), t, {ξ, ψ, 0}],
t > 0]
(*
E^(-(((x - ξ)^2 + (y - ψ)^2)/(4 d t)))/(4 d π t)
*)