程序1：用FDTD方法分析二維TM波在正弦波源激勵情況下，正方形Mur吸收邊界的性能
clear;clc;
N=100;
c=3e8;
f=1.5*10^(9);
lambda=c/f
eb=8.85e-12;
u=4*pi*10^-7;
Ez1=zeros(N+1,N+1);
Hx1=zeros(N+1,N);
Hy1=zeros(N,N+1);
Ez2=zeros(N+1,N+1);
Hx2=zeros(N+1,N);
Hy2=zeros(N,N+1);
deltaL=lambda/10;
deltaT=deltaL/(2*c);
 
for  n=1:200 
    ca=1;
    cb=deltaT/eb;
    cp=1;
    cq=deltaT/u;%迭代系數(shù)
    It=sin(2*pi*f*n*deltaT);
    %非邊界場值迭代
    clear i j;
    for i=1:(N+1)
        for j=1:N
            Hx2(i,j)=cp*Hx1(i,j)-cq*(Ez1(i,j+1)-Ez1(i,j))/deltaL;
        end
    end
    clear i j;
    for i=1:N
        for j=1:(N+1)
            Hy2(i,j)=cp*Hy1(i,j)+cq*(Ez1(i+1,j)-Ez1(i,j))/deltaL;
        end
    end
    clear i j;
    for i=2:N
        for j=2:N
Ez2(i,j)=ca*Ez1(i,j)+cb*((Hy2(i,j)-Hy2(i-1,j))/deltaL-(Hx2(i,j)-Hx2(i,j-1))/deltaL);
        end
    end
    clear i j;
    i=50;j=50;
    Ez2(i,j)=Ez1(i,j)+It;%激勵加入
    %邊界非角點場值迭代
%左邊界
    clear i j;
    i=1;
    for j=2:N
Ez2(i,j)=Ez1(i+1,j)+(c*deltaT-deltaL)*(Ez2(i+1,j)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%右邊界
    clear i j;
    i=N+1;
    for j=2:N
Ez2(i,j)=Ez1(i-1,j)+(c*deltaT-deltaL)*(Ez2(i-1,j)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%上邊界
    clear i j;
    j=N+1;
    for i=2:N
Ez2(i,j)=Ez1(i,j-1)+(c*deltaT-deltaL)*(Ez2(i,j-1)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%下邊界
    clear i j;
    j=1;
    for i=2:N
Ez2(i,j)=Ez1(i,j+1)+(c*deltaT-deltaL)*(Ez2(i,j+1)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%角點的迭代
    clear i j;
Ez2(1,1)=Ez1(2,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,2)-Ez1(1,1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(N+1,N+1)=Ez1(N,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,N)-Ez1(N+1,N+1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(1,N+1)=Ez1(2,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,N)-Ez1(1,N+1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(N+1,1)=Ez1(N,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,2)-Ez1(N+1,1))/(c*deltaT+sqrt(2)*deltaL);
%一次迭代結(jié)束，新得到的值變?yōu)橄麓蔚呐f值
    clear i j;
    for i=1:(N+1);
        for j=1:(N+1)
            Ez1(i,j)=Ez2(i,j);
        end
    end
    clear i j;
    for i=1:(N+1)
        for j=1:N
            Hx1(i,j)=Hx2(i,j);
        end
    end
     clear i j;
    for i=1:N
        for j=1:(N+1)
            Hy1(i,j)=Hy2(i,j);
        end
    end
%一次迭代結(jié)束
end
clear i j;
i=1:(N+1);
j=1:(N+1);
figure(1);
contour(i,j,Ez2(i,j));xlabel('x');ylabel('y');