matlab - Mathematica または Matlab による拡散方程式

Question

ご不便おかけしてすみません、

この特定の拡散方程式を NO 境界で解こうとしますhttp://astronomy.nju.edu.cn/~chenpf/c/courses/fluid/pringle81.pdf式 (2.10) nu=cost を使用します。このコードを使用して方程式を単純化します。

pde = D[s[x, t], 
    t] == (3/2)*D[s[x, t], x] + (3/4)*x^-2*D[s[x, t], {x, 2}] - 
    1/(4*x) ;
mu = 0.5;
sigma = 0.05;

そして方程式を解くために（nu = costを選択することによって、それは線形拡散部分方程式です）私はこれを使用します：

 sol = DSolve[{pde, 
   s[x, 0] == Exp[-(x - mu)^2/(2*sigma^2)]/Sqrt[2*Pi*sigma^2]}, 
  s[x, t], {x, t}]

初期関数（ガウス関数）の特定の選択をします。

しかし、私がそれをプロットしようとすると:

Plot3D[s[x, t] /. sol, {x, 0, 1}, {y, 0, Automatic}]

上記の論文の図 1 のプロットを再現するには、多くのエラーがあり、その理由がわかりません。

さらに、このMatlabコードを見つけました。このMatlabコードは、うまく機能するNO境界のある拡散タイプの方程式を再現しますが、方程式自体を変更して方程式のものを再現する方法を理解できません。上記の論文の (2.10)。

numx = 101;   %number of grid points in x
numt = 2000;  %number of time steps to be iterated
dx = 1/(numx - 1);
dt = 0.00005;

x = 0:dx:1;   %vector of x values, to be used for plotting

C = zeros(numx,numt);   %initialize everything to zero

%specify initial conditions
t(1) = 0;      %t=0
mu = 0.5;      
sigma = 0.05;
for i=1:numx
   C(i,1) = exp(-(x(i)-mu)^2/(2*sigma^2)) / sqrt(2*pi*sigma^2);
end

%iterate difference equations
for j=1:numt
   t(j+1) = t(j) + dt;
   for i=2:numx-1
      C(i,j+1) = C(i,j) + (dt/dx^2)*(C(i+1,j) - 2*C(i,j) + C(i-1,j)); 
   end
   C(1,j+1) = C(2,j+1);          %C(1,j+1) found from no-flux condition
   C(numx,j+1) = C(numx-1,j+1);  %C(numx,j+1) found from no-flux condition
end

figure(1);
hold on;
plot(x,C(:,1));
plot(x,C(:,11));
plot(x,C(:,101));
plot(x,C(:,1001));
plot(x,C(:,2001));
xlabel('x');
ylabel('c(x,t)');

誰か助けてくれませんか？

どうもありがとう。