利用时域有限差分FDTD仿真三维3D电磁波的传播，吸收边界条件为完全匹配层PML。

%***********************************************************************
%     3-D FDTD code with PML absorbing boundary conditions
%***********************************************************************
%                      Lijijun 
%                      Department of Physics  Ynagtze University
%                      lijijun@163.net
%                   
%    Date of this version: August 2006                  
%
%***********************************************************************
function fdtd3D

clear
clc
%***********************************************************************
%     Fundamental constants
%***********************************************************************
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/（cc*cc*muz）;       %permittivity of free space
etaz=sqrt（muz/epsz）;

freq=1e12;
lambda=cc/freq;
omega=2.0*pi*freq;    

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=101;           %number of grid cells in x-direction
if mod（ie2）==0
   ie=ie+1;
end       %Force it to be a odd integer!
je=131;            %number of grid cells in y-direction
if mod（je2）==0
   je=je+1;
end       %Force it to be a odd integer!
ke=51;             %number of grid cells in z-direction
if mod（ke2）==0
   ke=ke+1;
end       %Force it to be a odd integer!

ib=ie+1;
jb=je+1;
kb=ke+1;

is=（ie+1）/2;            %location of  hard source in x axis
js=（je+1）/2;          %location of  hard source in y axis
ks=（ke+1）/2;          %location of  hard source in z axis

ds=lambda/20;        %space increment of square lattice
dt=ds/（2.0*cc）;   %time step

nmax=100;%         %total number of time steps

iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
kebc=8;           %thickness of bottom and top PML region
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
kefbc=ke+2*kebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
kbfbc=kefbc+1;

rmax=1.0e-7;
orderbc=2;

%***********************************************************************
%     Material parameters
%***********************************************************************

media=2;

eps=[1.0 1.0];
sig=[0.0 1.0e+7];
mur=[1.0 1.0];
sim=[0.0 0.0];

ds = ds/sqrt（eps（1）*mur（1））;

%***********************************************************************
%     Wave excitation
%***********************************************************************

source=zeros（1nmax）;

for n=1:nmax
    source（n）=sin（omega*n*dt）;
end

%***********************************************************************
%     Field arrays
%***********************************************************************

ex=zeros（iejbkb）;           %fields in main grid 
ey=zero