% Physical constantseps0 = 8.8541878e-12;         % Permittivity of vacuummu0  = 4e-7 * pi;             % Permeability of vacuumc0   = 299792458;             % Speed of light in vacuum% Parameter initiationLx = .05; Ly = .04; Lz = .03; % Cavity dimensions in metersNx =  25; Ny =  20; Nz =  15; % Number of cells in each directionCx = Nx / Lx;                 % Inverse cell dimensionsCy = Ny / Ly;Cz = Nz / Lz;Nt = 32;                    % Number of time stepsDt = 1/(c0*norm([Cx Cy Cz])); % Time step% Allocate field matricesEx = zeros(Nx  , Ny+1, Nz+1);Ey = zeros(Nx+1, Ny  , Nz+1);Ez = zeros(Nx+1, Ny+1, Nz  );Hx = zeros(Nx+1, Ny  , Nz  );Hy = zeros(Nx  , Ny+1, Nz  );Hz = zeros(Nx  , Ny  , Nz+1);% Allocate time signalsEt = zeros(Nt,3);% Initiate fields (near but not on the boundary)Ex(1,2,2) = 1;Ey(2,1,2) = 2;Ez(2,2,1) = 3;% Time steppingfor n = 1:Nt;  % Update H everywhere   for i = 1:Nx+1     for j = 1:Ny        for k = 1:Nz           Hx(i,j,k) = Hx(i,j,k) + (Dt/mu0)* ...           ((Ey(i,j,k+1)-Ey(i,j,k))*Cz - (Ez(i,j+1,k)-Ez(i,j,k))*Cy);        end     end  end    Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,:)-Ez(1:Nx,:,:))*Cx ...                    - (Ex(:,:,2:Nz+1)-Ex(:,:,1:Nz))*Cz);  Hz = Hz + (Dt/mu0)*((Ex(:,2:Ny+1,:)-Ex(:,1:Ny,:))*Cy ...                    - (Ey(2:Nx+1,:,:)-Ey(1:Nx,:,:))*Cx);% Update E everywhere except on boundary  for i = 1:Nx     for j = 2:Ny        for k = 2:Nz           Ex(i,j,k) = Ex(i,j,k) + (Dt /eps0) * ...              ((Hz(i,j,k)-Hz(i,j-1,k))*Cy-(Hy(i,j,k)-Hy(i,j,k-1))*Cz);        end     end  end  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...      ((Hx(2:Nx,:,2:Nz)-Hx(2:Nx,:,1:Nz-1))*Cz ...     - (Hz(2:Nx,:,2:Nz)-Hz(1:Nx-1,:,2:Nz))*Cx);  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...      ((Hy(2:Nx,2:Ny,:)-Hy(1:Nx-1,2:Ny,:))*Cx ...     - (Hx(2:Nx,2:Ny,:)-Hx(2:Nx,1:Ny-1,:))*Cy);  % Sample the electric field at chosen points  Etss(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];end