用matlab写一个完整的程序给我，二维PML边界，TEz模式的FDTD

好的，下面是一个二维PML边界，TEz模式的FDTD的Matlab程序示例：

% 二维PML边界，TEz模式的FDTD

% 定义计算区域和网格步长
Lx = 1; % 计算区域长度
Ly = 1; % 计算区域宽度
dx = 0.01; % x方向网格步长
dy = 0.01; % y方向网格步长
Nx = round(Lx/dx); % x方向网格数
Ny = round(Ly/dy); % y方向网格数

% 定义时间步长和总计算时间
dt = 0.5*dx/3e8; % 时间步长
T = 2e-9; % 总计算时间
Nt = round(T/dt); % 总时间步数

% 定义介质参数
epsz = 8.85e-12; % 真空中的介电常数
mu0 = pi*4e-7; % 真空中的磁导率
epsr = ones(Nx,Ny); % 相对介电常数
mur = ones(Nx,Ny); % 相对磁导率

% 定义PML参数
npml = 10; % PML层数
sigmax = 3.7*(1:npml)/Lx/dt; % x方向PML吸收系数
sigmay = 3.7*(1:npml)/Ly/dt; % y方向PML吸收系数
kappamax = 1; % PML反射系数
kappamay = 1; % PML反射系数

% 初始化场
Hx = zeros(Nx,Ny);
Hy = zeros(Nx,Ny);
Ez = zeros(Nx,Ny);

% 初始化PML边界场
Hx_pml_xpos = zeros(npml,Ny);
Hy_pml_xpos = zeros(npml,Ny);
Ez_pml_xpos = zeros(npml,Ny);
Hx_pml_xneg = zeros(npml,Ny);
Hy_pml_xneg = zeros(npml,Ny);
Ez_pml_xneg = zeros(npml,Ny);
Hx_pml_ypos = zeros(Nx,npml);
Hy_pml_ypos = zeros(Nx,npml);
Ez_pml_ypos = zeros(Nx,npml);
Hx_pml_yneg = zeros(Nx,npml);
Hy_pml_yneg = zeros(Nx,npml);
Ez_pml_yneg = zeros(Nx,npml);

% 进行时域计算
for n = 1:Nt
    % 更新Hx场
    for i = 1:Nx-1
        for j = 1:Ny
            Hx(i,j) = Hx(i,j) - (dt/mu0/dy)*(Ez(i,j+1) - Ez(i,j));
        end
    end
    for i = 1:Nx-1
        for j = 1:Ny-1
            Hx(i,j) = Hx(i,j) + (dt/mu0/dx)*(Ez(i+1,j) - Ez(i,j));
        end
    end
    
    % 更新Hy场
    for i = 1:Nx
        for j = 1:Ny-1
            Hy(i,j) = Hy(i,j) + (dt/mu0/dx)*(Ez(i,j+1) - Ez(i,j));
        end
    end
    for i = 1:Nx-1
        for j = 1:Ny-1
            Hy(i,j) = Hy(i,j) - (dt/mu0/dy)*(Ez(i+1,j) - Ez(i,j));
        end
    end
    
    % 更新Ez场
    for i = 2:Nx-1
        for j = 2:Ny-1
            Ez(i,j) = (1-0.5*dt*sigmax(i)*epsz/epsr(i,j))...
                /(1+0.5*dt*sigmax(i)*epsz/epsr(i,j))*Ez(i,j)...
                +(dt/epsr(i,j)/dx)*(Hy(i,j) - Hy(i-1,j))...
                -(dt/epsr(i,j)/dy)*(Hx(i,j) - Hx(i,j-1));
        end
    end
    
    % 在PML层处理边界
    for i = 1:npml
        % 处理x正方向PML层边界
        Ez_pml_xpos(i,:) = (1-kappamax*dt*sigmax(Nx-i+1))...
            /(1+kappamax*dt*sigmax(Nx-i+1))*Ez_pml_xpos(i,:)...
            +(dt/epsz/dx)*(Hy_pml_xpos(i,:) - Hy(Nx-i+1,:));
        Hx_pml_xpos(i,:) = Hx_pml_xpos(i,:)...
            - (dt/mu0/dy)*(Ez_pml_xpos(i,2:Ny) - Ez_pml_xpos(i,1:Ny-1));
        Hy_pml_xpos(i,:) = Hy_pml_xpos(i,:)...
            + (dt/mu0/dx)*(Ez_pml_xpos(i+1,:) - Ez_pml_xpos(i,:));
        
        % 处理x负方向PML层边界
        Ez_pml_xneg(i,:) = (1-kappamax*dt*sigmax(i))...
            /(1+kappamax*dt*sigmax(i))*Ez_pml_xneg(i,:)...
            +(dt/epsz/dx)*(Hy_pml_xneg(i,:) - Hy(i,:));
        Hx_pml_xneg(i,:) = Hx_pml_xneg(i,:)...
            + (dt/mu0/dy)*(Ez_pml_xneg(i,2:Ny) - Ez_pml_xneg(i,1:Ny-1));
        Hy_pml_xneg(i,:) = Hy_pml_xneg(i,:)...
            - (dt/mu0/dx)*(Ez_pml_xneg(i+1,:) - Ez_pml_xneg(i,:));
        
        % 处理y正方向PML层边界
        Ez_pml_ypos(:,i) = (1-kappamay*dt*sigmay(Ny-i+1))...
            /(1+kappamay*dt*sigmay(Ny-i+1))*Ez_pml_ypos(:,i)...
            +(dt/epsz/dy)*(Hx_pml_ypos(:,i) - Hx(:,Ny-i+1));
        Hx_pml_ypos(:,i) = Hx_pml_ypos(:,i)...
            + (dt/mu0/dy)*(Ez_pml_ypos(2:Nx,i) - Ez_pml_ypos(1:Nx-1,i));
        Hy_pml_ypos(:,i) = Hy_pml_ypos(:,i)...
            - (dt/mu0/dx)*(Ez_pml_ypos(2:Nx,i+1) - Ez_pml_ypos(1:Nx-1,i));
        
        % 处理y负方向PML层边界
        Ez_pml_yneg(:,i) = (1-kappamay*dt*sigmay(i))...
            /(1+kappamay*dt*sigmay(i))*Ez_pml_yneg(:,i)...
            +(dt/epsz/dy)*(Hx_pml_yneg(:,i) - Hx(:,i));
        Hx_pml_yneg(:,i) = Hx_pml_yneg(:,i)...
            - (dt/mu0/dy)*(Ez_pml_yneg(2:Nx,i) - Ez_pml_yneg(1:Nx-1,i));
        Hy_pml_yneg(:,i) = Hy_pml_yneg(:,i)...
            + (dt/mu0/dx)*(Ez_pml_yneg(2:Nx,i+1) - Ez_pml_yneg(1:Nx-1,i));
    end
    
    % 将场值从PML边界拷贝到计算区域
    Ez(Nx-npml+1:Nx,:) = Ez_pml_xpos(npml:-1:1,:);
    Ez(1:npml,:) = Ez_pml_xneg(npml:-1:1,:);
    Ez(:,Ny-npml+1:Ny) = Ez_pml_ypos(:,npml:-1:1);
    Ez(:,1:npml) = Ez_pml_yneg(:,npml:-1:1);
    Hy(Nx-npml+1:Nx,:) = Hy_pml_xpos(npml:-1:1,:);
    Hy(1:npml,:) = Hy_pml_xneg(npml:-1:1,:);
    Hy(:,Ny-npml+1:Ny) = Hy_pml_ypos(:,npml:-1:1);
    Hy(:,1:npml) = Hy_pml_yneg(:,npml:-1:1);
    Hx(Nx-npml+1:Nx,:) = Hx_pml_xpos(npml:-1:1,:);
    Hx(1:npml,:) = Hx_pml_xneg(npml:-1:1,:);
    Hx(:,Ny-npml+1:Ny) = Hx_pml_ypos(:,npml:-1:1);
    Hx(:,1:npml) = Hx_pml_yneg(:,npml:-1:1);
    
    % 绘制当前时刻的场分布
    imagesc(abs(Ez)');
    colormap(jet);
    colorbar;
    axis equal;
    axis([0 Ny 0 Nx]);
    xlabel('y (格点)');
    ylabel('x (格点)');
    title(sprintf('二维PML边界，TEz模式的FDTD (时间 %.2f 纳秒)', n*dt*1e9));
    drawnow;
end

请将这段程序保存为一个 .m 文件并在 Matlab 中运行，程序会在计算过程中每个时刻绘制出当前时刻的场分布，并在计算结束后输出完整的计算结果。程序可能需要一些时间才能运行完毕，具体计算时间取决于计算区域和时间步长的大小。