程序代碼如下：
% 二維FDTD  TE波仿真 
clear all;
% 定義常數(shù)
pi=3.1415;
c=3.0e10;                     %高斯制下光速
f=1.0e15;                     %頻率
lambda=c/f;                   %波長
nmax=300;                     %時(shí)間步數(shù)
del_s=lambda/20;              %每最小波長20個(gè)采樣點(diǎn)          
del_t=0.5*del_s/c;            %迭代時(shí)間步長
n=182;                        %真空區(qū)域網(wǎng)格數(shù)
np=8;                         %pml層數(shù)
N1=n+2*np;                    %總網(wǎng)格數(shù)
N=N1+1;                       %采樣點(diǎn)數(shù)
M=4;                          %導(dǎo)電率漸變指數(shù)
sigma_max=(M+1)/1.50/pi/del_s;%最大導(dǎo)電率

% TE波的分量初始化
tic;
figure(1);
axis([0 N 0 N -0.5 0.5]);
Ex=zeros(N1,N);               %x方向?yàn)闄M向，采樣點(diǎn)為網(wǎng)格的橫向邊，故行數(shù)+1
Ey=zeros(N,N1);               %y方向?yàn)榭v向，采樣點(diǎn)為網(wǎng)格的縱向邊，故列數(shù)+1
Bz=zeros(N1,N1);              %矩陣行為縱向網(wǎng)格數(shù)，矩陣列為橫向網(wǎng)格數(shù),循環(huán)中用j表示行數(shù)，i表示列數(shù)
Bzx=zeros(N1,N1);
Bzy=zeros(N1,N1);
Bzxx=zeros(nmax,2);

%進(jìn)入電磁場迭代計(jì)算
for tt=1:nmax
  for i=1:N1
        if i>=np+1&&i<=N1-np
           di=0;
        elseif i<=np
           di=np-i+0.5;
        elseif i>=N1-np+1
           di=np+i-N1-0.5;
        end                                 %di是采樣點(diǎn)橫向距PML內(nèi)邊界的距離
         sigma_mx=sigma_max*(di/np)^M;
     for j=1:N1  
         if j>=np+1&&j<=N1-np
                dj=0;
         elseif j<=np
                dj=np-j+0.5;
         elseif j>=N1-np+1
                dj=np+j-N1-0.5;
         end                                 %dj是采樣點(diǎn)縱向距PML內(nèi)邊界的距離
            sigma_my=sigma_max*(dj/np)^M;                                     
         
              if sigma_mx+sigma_my==0        %真空區(qū)
                  if j==100&&i==100
                      t=30;                %可選擇高斯脈沖
                      term=(tt-t);
                      pulse=exp(-4*pi*term^2/20^2) ;
%                       pulse=sin(2*pi*tt/40); %可選正弦時(shí)諧源
                      c_miu=c*del_t/del_s;
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脈沖源
                  else
                      c_miu=c*del_t/del_s;
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2;
                  end
              else                            %PML區(qū)
                    cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t);
                    cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s;
                    Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));               
                    cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t);
                    cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s;
                    Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j));
                    Bz(i,j)=Bzx(i,j)+Bzy(i,j);
              end
     end
  end
  for i=2:N1
      if i>=np+1&&i<=N1-np
         di=0;
      elseif i<=np
         di=np-i+1;
      elseif i>=N1-np+1
         di=np+i-N1-1;
      end                                   %di是采樣點(diǎn)橫向距PML內(nèi)邊界的距離
      sigma_ex=sigma_max*(di/np)^M;
     for j=1:N1
         cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);
         cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s;
         Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j));
     end
  end
  for i=1:N1
      for j=2:N1 
          if j>=np+1&&j<=N1-np
             dj=0;
          elseif j<=np
             dj=np-j+1;
          elseif j>=N1-np+1
             dj=np+j-N1-1;
          end                               %dj是采樣點(diǎn)縱向距PML內(nèi)邊界的距離
          sigma_ey=sigma_max*(dj/np)^M;
          cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t);
          cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;
          Ex(i,j)=cam*Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1));   
      end 
  end    
  Bzxx(tt,1)=Bz(90,50);                     %對靠近邊界的中央磁場點(diǎn)采樣
  Bzxx(tt,2)=Bz(50,90);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%電磁場的計(jì)算部分結(jié)束
  figure(1);                                %可視化處理
  clf;   
  mesh(Bz);                                 %磁場的幅值
  axis([0 N 0 N -0.5 0.5]);
  xlabel('i')
  ylabel('j')
  drawnow;
end
figure(2);
plot(abs(Bzxx));
figure(3);
title('磁場幅值分布圖');
surface(abs(Bz).^2);
shading interp;
axis square;
toc