% 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 = 1024;                    % 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  Hx = Hx + (Dt/mu0)*(diff(Ey,1,3)*Cz - diff(Ez,1,2)*Cy);  Hy = Hy + (Dt/mu0)*(diff(Ez,1,1)*Cx - diff(Ex,1,3)*Cz);  Hz = Hz + (Dt/mu0)*(diff(Ex,1,2)*Cy - diff(Ey,1,1)*Cx);  % Update E everywhere except on boundary  Ex(:,2:Ny,2:Nz) = Ex(:,2:Ny,2:Nz) + (Dt /eps0) * ...      (diff(Hz(:,:,2:Nz),1,2)*Cy - diff(Hy(:,2:Ny,:),1,3)*Cz);  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...      (diff(Hx(2:Nx,:,:),1,3)*Cz - diff(Hz(:,:,2:Nz),1,1)*Cx);  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...      (diff(Hy(:,2:Ny,:),1,1)*Cx - diff(Hx(2:Nx,:,:),1,2)*Cy);  % Sample the electric field at chosen points%   Eto(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];   Eto(n,:) = [Ex(10,10,10) Ey(10,10,10) Ez(10,10,10)];  end%===============================================================Eto1(1:1024)=sqrt(Eto(1:1024,1).^2+Eto(1:1024,2).^2+Eto(1:1024,3).^2);% plot(Eto(:,1),':');% hold on% plot(Eto(:,2),'-.');% plot(Eto(:,3),'--');% plot(Eto1,'o');plot(Eto(:,1));hold onplot(Eto(:,2),'m');plot(Eto(:,3),'g');plot(Eto1,'r');