My 3D FDTD (UPML) MATLAB code is not updating and I need help.

[zozokim]

I've been studying the UPML portion of Computational electrodynamics by ALLEN TAFLOVE and writing matlab code, but the electromagnetic field is not updating.
Please help me

I can't remove things like emoticon :/), so I'm attaching it as a notepad.
Thank you.

I apologize for my poor English.
This is my MATLAB code.

--------------------------------------------------------------------------------------------------------

clear all; close all; clc;

%%
%Defining constants

e0 = 8.854*10^-12;                      % Permittivity of vacuum [F/m]
u0 = 4*pi*10^-7;                        % Permeability of vacuum [H/m]
c0 = 1/(e0*u0)^0.5;                     % Speed of light         [m/s]
f0 = 10e6;                              % Frequency of Source    [Hz]
w0 = 2*pi*f0;                           % Angular frequency      [rad/s]
l0 = c0/f0;                             % lambda                 [m]

X = 100;                                % X size [m]
Y = 100;                                % Y size [m]
Z = 100;                                % Z size [m]
Mx = 100;                               % X cell
My = 100;                               % Y cell
Mz = 100;                               % Z cell
N = 2000;                                % Time step
dx = X/Mx;                              % Increment of X
dy = Y/My;                              % Increment of Y
dz = Z/Mz;                              % Increment of Z
dt = (dx^-2+dy^-2+dz^-2)^(-0.5)/c0;    % Courant number (Magic time step)
x_axis = (1:Mx)*dx;                     % Physical length of X
y_axis = (1:My)*dy;                     % Physical length of X


sigx = zeros(Mx+1,My+1,Mz+1);                 % Conductivity x-axis
sigy = zeros(Mx+1,My+1,Mz+1);                 % Conductivity y-axis
sigz = zeros(Mx+1,My+1,Mz+1);                 % Conductivity z-axis (In 2D, sigz = 0)
e = e0*ones(Mx+1,My+1,Mz+1);                  % Permittivity matrix of PML region
u = u0*ones(Mx+1,My+1,Mz+1);                  % Permeability matrix of PML region
eta = sqrt(u0/e0);                      % Intrinsic impedance of PML region

%PML condition(polynomial grading)

p = 25;                                 % PML grid
kx = 1;                                 % Kappa x
ky = 1;                                 % Kappa y
kz = 1;                                 % Kappa z (In 2D, Kappa z = 1)
m = 4;                                  % Grading order
R0 = exp(-16);                          % Reflection errors
sig_max = -((m+1)*log(R0))/(2*eta*p);   % Sigma max

for n = 1:p
    sigx(p+1-n,:,:) = (n/p)^m*sig_max;
    sigy(:,p+1-n,:) = (n/p)^m*sig_max;
    sigz(:,:,p+1-n) = (n/p)^m*sig_max;
    sigx(end-p+n,:,:) = (n/p)^m*sig_max;
    sigy(:,end-p+n,:) = (n/p)^m*sig_max;
    sigz(:,:,end-p+n) = (n/p)^m*sig_max;
end

%Initializing dimensions
Dx = zeros(Mx,My+1,Mz+1);
Ex = zeros(Mx,My+1,Mz+1);
Dy = zeros(Mx+1,My,Mz+1);
Ey = zeros(Mx+1,My,Mz+1);
Dz = zeros(Mx+1,My+1,Mz);
Ez = zeros(Mx+1,My+1,Mz);
Bx = zeros(Mx+1,My,Mz);
Hx = zeros(Mx+1,My,Mz);
By = zeros(Mx,My+1,Mz);
Hy = zeros(Mx,My+1,Mz);
Bz = zeros(Mx,My,Mz+1);
Hz = zeros(Mx,My,Mz+1);

% monitor = zeros(Mx,My,N);

Wx = 2*e0*kx-sigx*dt;                   %Coefficient of update equations
Wy = 2*e0*ky-sigy*dt;                   %Coefficient of update equations
Wz = 2*e0*kz-sigz*dt;                   %Coefficient of update equations
Qx = 2*e0*kx+sigx*dt;                   %Coefficient of update equations
Qy = 2*e0*ky+sigy*dt;                   %Coefficient of update equations
Qz = 2*e0*kz+sigz*dt;                   %Coefficient of update equations

clims = [-70 -30];

% v = VideoWriter('2D FDTD PML1');
% open(v);

%%
%UPML Update equation by auxiliary differential equation

% Source
t = (1:N)*dt;
E_in = exp(-2 * pi^2 * (f0)^2 * (t - 1/(f0)).^2 );
figure(2); 
plot(t, E_in);

i1 = 1:Mx;
j1 = 1:My;
k1 = 1:Mz;
i2 = 2:Mx;
j2 = 2:My;
k2 = 2:Mz;

% Define coefficients of update equations
T1 = (2*e0*dt);

Dx_C1 = Wy(i1,j2,k2)./Qy(i1,j2,k2);
Dx_C2 = T1./Qy(i1,j2,k2)/dy;
Dx_C3 = T1./Qy(i1,j2,k2)/dz;
Dy_C1 = Wz(i2,j1,k2)./Qz(i2,j1,k2);
Dy_C2 = T1./Qz(i2,j1,k2)/dz;
Dy_C3 = T1./Qz(i2,j1,k2)/dx;
Dz_C1 = Wx(i2,j2,k1)./Qx(i2,j2,k1);
Dz_C2 = T1./Qx(i2,j2,k1)/dx;
Dz_C3 = T1./Qx(i2,j2,k1)/dy;

Ex_C1 = Wz(i1,j2,k2)./Qz(i1,j2,k2);
Ex_C2 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Qx(i1,j2,k2);
Ex_C3 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Wx(i1,j2,k2);
Ey_C1 = Wx(i2,j1,k2)./Qx(i2,j1,k2);
Ey_C2 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Qy(i2,j1,k2);
Ey_C3 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Wy(i2,j1,k2);
Ez_C1 = Wy(i2,j2,k1)./Qy(i2,j2,k1);
Ez_C2 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Qz(i2,j2,k1);
Ez_C3 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Wz(i2,j2,k1);

Bx_C1 = Wy(:,j1,k1)./Qy(:,j1,k1);
Bx_C2 = T1./Qy(:,j1,k1)/dy;
Bx_C3 = T1./Qy(:,j1,k1)/dz;
By_C1 = Wz(i1,:,k1)./Qz(i1,:,k1);
By_C2 = T1./Qz(i1,:,k1)/dz;
By_C3 = T1./Qz(i1,:,k1)/dx;
Bz_C1 = Wx(i1,j1,:)./Qx(i1,j1,:);
Bz_C2 = T1./Qx(i1,j1,:)/dx;
Bz_C3 = T1./Qx(i1,j1,:)/dy;

Hx_C1 = Wz(:,j1,k1)./Qz(:,j1,k1);
Hx_C2 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Qx(:,j1,k1);
Hx_C3 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Wx(:,j1,k1);
Hy_C1 = Wx(i1,:,k1)./Qx(i1,:,k1);
Hy_C2 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Qy(i1,:,k1);
Hy_C3 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Wy(i1,:,k1);
Hz_C1 = Wy(i1,j1,:)./Qy(i1,j1,:);
Hz_C2 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Qz(i1,j1,:);
Hz_C3 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Wz(i1,j1,:);


figure(1);
for n = 1:N
    t(n)
    
    Dx_temp = Dx(:,j2,k2);
    Dy_temp = Dy(i2,:,k2);
    Dz_temp = Dz(i2,j2,:);
    
    Ex_temp = Ex(:,j2,k2);
    Ey_temp = Ey(i2,:,k2);
    Ez_temp = Ez(i2,j2,:);
    
    Bx_temp = Bx(:,j1,k1);
    By_temp = By(i1,:,k1);
    Bz_temp = Bz(i1,j1,:);
    
    Hx_temp = Hx(:,j1,k1);
    Hy_temp = Hy(i1,:,k1);
    Hz_temp = Hz(i1,j1,:);
    
    Dx(:,j2,k2) = Dx_C1.*Dx_temp + Dx_C2.*(Hz(:,j2,k2) - Hz(:,j2-1,k2)) - Dx_C3.*(Hy(:,j2,k2) - Hy(:,j2,k2-1));
    Dy(i2,:,k2) = Dy_C1.*Dy_temp + Dy_C2.*(Hx(i2,:,k2) - Hx(i2,:,k2-1)) - Dy_C3.*(Hz(i2,:,k2) - Hz(i2-1,:,k2));
    Dz(i2,j2,:) = Dz_C1.*Dz_temp + Dz_C2.*(Hy(i2,j2,:) - Hy(i2-1,j2,:)) - Dz_C3.*(Hx(i2,j2,:) - Hx(i2,j2-1,:));
    
    Ex(:,j2,k2) = Ex_C1.*Ex_temp + Ex_C2.*Dx(:,j2,k2) - Ex_C3.*Dx_temp;
    Ey(i2,:,k2) = Ey_C1.*Ey_temp + Ey_C2.*Dy(i2,:,k2) - Ey_C3.*Dy_temp;
    Ez(i2,j2,:) = Ez_C1.*Ez_temp + Ez_C2.*Dz(i2,j2,:) - Ez_C3.*Dz_temp;
    
    Bx(:,j1,k1) = Bx_C1.*Bx_temp - Bx_C2.*(Ez(:,j1+1,k1) - Ez(:,j1,k1)) + Bx_C3.*(Ey(:,j1,k1+1) - Ey(:,j1,k1));
    By(i1,:,k1) = By_C1.*By_temp - By_C2.*(Ex(i1,:,k1+1) - Ex(i1,:,k1)) + By_C3.*(Ez(i1+1,:,k1) - Ez(i1,:,k1));
    Bz(i1,j1,:) = Bz_C1.*Bz_temp - Bz_C2.*(Ey(i1+1,j1,:) - Ey(i1,j1,:)) + Bz_C3.*(Ex(i1,j1+1,:) - Ex(i1,j1,:));
    
    Hx(:,j1,k1) = Hx_C1.*Hx_temp + Hx_C2.*Bx(:,j1,k1) - Hx_C3.*Bx_temp;
    Hy(i1,:,k1) = Hy_C1.*Hy_temp + Hy_C2.*By(i1,:,k1) - Hy_C3.*By_temp;
    Hz(i1,j1,:) = Hz_C1.*Hz_temp + Hz_C2.*Bz(i1,j1,:) - Hz_C3.*Bz_temp;
    
    Ez(round(X/2),round(Y/2),round(Z/2)) = Ez(round(X/2),round(Y/2),round(Z/2)) - exp(-2 * pi^2 * (f0)^2 * (n*dt - 1/(f0)).^2 );
    
    E_total = sqrt(abs(Ex(i1,j1,k1)).^2+abs(Ey(i1,j1,k1)).^2+abs(Ez(i1,j1,k1)).^2);
%     subplot(1,2,1)
    %pcolor(db(E_total(:,:,round(Mz/2))))
    figure(1);
    pcolor(y_axis, x_axis, (E_total(:,:,round(Z/2))))
    shading flat
    caxis([0 5]);
    colorbar;
    
    drawnow;
%     frame = getframe(gcf);
%     writeVideo(v,frame);
end

% close(v);

---------------------------------------------------------------------------------------------------------------------------