音響FDTD1次元バージョン

MATLABで音響FDTDプログラムを作成しました．とりあえず1次元バージョンです．
*自由音場における音圧p，粒子速度uに関する偏微分方程式を差分化して，逐一音圧と粒子速度を交互に計算していきます．
*密度と体積弾性率は室温20度の値を用いています．
*Murの1次吸収境界を用いています．
*初期波源は音圧に1次元ガウシアンパルスを与え，粒子速度をゼロとしています．

%===Parameters===%
NT=200;NX=101;
SIG=10e-3/sqrt(2);
RHO=1.21;KKK=1.42e5;
cc=sqrt(KKK/RHO);
dx=1e-3;dt=1e-6;
mur=(cc*dt-dx)/(cc*dt+dx);
duu=dt/RHO/dx;
dpp=KKK*dt/dx;
x=0:dx:(NX-1)*dx;
p=exp(-((x-0.05).^2)/2/SIG/SIG);
u=zeros(size(x));
for T=0:NT;
    t=T*dt;
    %===update u===%
    u(1:NX-1)=u(1:NX-1)-duu*(p(2:NX)-p(1:NX-1));
    %===update p===%
    tmp1=p(2);tmp2=p(NX-1);
    p(2:NX-1)=p(2:NX-1)-dpp*(u(2:NX-1)-u(1:NX-2));
    %===Mur1st===%
    p(1)=tmp1-mur*(p(1)-p(2));
    p(NX)=tmp2-mur*(p(NX)-p(NX-1));
    %===Plot===%
    plot(x,p);
    axis([0,max(x),-0.2,1.2]);
    drawnow;
end