Hy(1,:)=expboundary*Hy(1,:)+(1-expboundary)*...
      (Ez(1,:)-EzxPMLC(NPML,:)-EzyPMLC(NPML,:))/(SigmaBound*factor*Dx);  %Boundary C
   Hy(NTx,:)=expboundary*Hy(NTx,:)+(1-expboundary)*...
      (EzxPMLD(1,:)+EzyPMLD(1,:)-Ez(NTx-1,:))/(SigmaBound*factor*Dx);     %Boundary D
   
   %The following part is for ZONE A. H components.
   HxPMLA(:,1:NPML-1)=exp(-Sigmay_xA(:,1:NPML-1)*Dt/e0).*HxPMLA(:,1:NPML-1)-...
      (1-exp(-Sigmay_xA(:,1:NPML-1)*Dt/e0))./(Sigmay_xA(:,1:NPML-1)*factor*Dy).*...
      (EzxPMLA(:,2:NPML)+EzyPMLA(:,2:NPML)-EzxPMLA(:,1:NPML-1)-EzyPMLA(:,1:NPML-1));
   
   HyPMLA(2:NTx-1,:)=HyPMLA(2:NTx-1,:)+...
      Dt*(EzxPMLA(2:NTx-1,:)+EzyPMLA(2:NTx-1,:)-EzxPMLA(1:NTx-2,:)-EzyPMLA(1:NTx-2,:))/(mu0*Dx);
   
   HyPMLA(1,:)=expboundary*HyPMLA(1,:)+(1-expboundary)*...
      (EzxPMLA(1,:)+EzyPMLA(1,:)-EzxPML1(NPML,:)-EzyPML1(NPML,:))/(SigmaBound*factor*Dx);  %Boundary Left
   HyPMLA(NTx,:)=expboundary*HyPMLA(NTx,:)+(1-expboundary)*...
      (EzxPML2(1,:)+EzyPML2(1,:)-EzxPMLA(NTx-1,:)-EzyPMLA(NTx-1,:))/(SigmaBound*factor*Dx);  %Boundary Right
   HxPMLA(:,NPML)=HxPMLA(:,NPML-1); %Boundary Upper;
   
   %The following part is for ZONE B. 
   %H components.
   HxPMLB(:,2:NPML)=exp(-Sigmay_xB(:,2:NPML)*Dt/e0).*HxPMLB(:,2:NPML)-...
      (1-exp(-Sigmay_xB(:,2:NPML)*Dt/e0))./(Sigmay_xB(:,2:NPML)*factor*Dy).*...
      (EzxPMLB(:,2:NPML)+EzyPMLB(:,2:NPML)-EzxPMLB(:,1:NPML-1)-EzyPMLB(:,1:NPML-1));
   
   HyPMLB(2:NTx-1,:)=HyPMLB(2:NTx-1,:)+...
      Dt*(EzxPMLB(2:NTx-1,:)+EzyPMLB(2:NTx-1,:)-EzxPMLB(1:NTx-2,:)-EzyPMLB(1:NTx-2,:))/(mu0*Dx);
   
   HyPMLB(1,:)=expboundary*HyPMLB(1,:)+(1-expboundary)*...
      (EzxPMLB(1,:)+EzyPMLB(1,:)-EzxPML3(NPML,:)-EzyPML3(NPML,:))/(SigmaBound*factor*Dx);  %Boundary Left
   HyPMLB(NTx,:)=expboundary*HyPMLB(NTx,:)+(1-expboundary)*...
      (EzxPML4(1,:)+EzyPML4(1,:)-EzxPMLB(NTx-1,:)-EzyPMLB(NTx-1,:))/(SigmaBound*factor*Dx);  %Boundary Right
   HxPMLB(:,1)=HxPMLB(:,2); %Boundary Bottom;

   %The following part is for ZONE C. 
   %H components.
   HxPMLC(:,2:NTy-1)=HxPMLC(:,2:NTy-1)-...
      Dt*(EzxPMLC(:,2:NTy-1)+EzyPMLC(:,2:NTy-1)-EzxPMLC(:,1:NTy-2)-EzyPMLC(:,1:NTy-2))/(mu0*Dy);
   
   HyPMLC(2:NPML,:)=exp(-Sigmax_yC(2:NPML,:)*Dt/e0).*HyPMLC(2:NPML,:)+...
      (1-exp(-Sigmax_yC(2:NPML,:)*Dt/e0))./(Sigmax_yC(2:NPML,:)*factor*Dx).*...
      (EzxPMLC(2:NPML,:)+EzyPMLC(2:NPML,:)-EzxPMLC(1:NPML-1,:)-EzyPMLC(1:NPML-1,:));
   
   HxPMLC(:,1)=expboundary*HxPMLC(:,1)-(1-expboundary)*...
      (EzxPMLC(:,1)+EzyPMLC(:,1)-EzxPML3(:,NPML)-EzyPML3(:,NPML))/(SigmaBound*factor*Dy);  %Boundary Down
   HxPMLC(:,NTy)=expboundary*HxPMLC(:,NTy)-(1-expboundary)*...
      (EzxPML1(:,1)+EzyPML1(:,1)-EzxPMLC(:,NTy-1)-EzyPMLC(:,NTy-1))/(SigmaBound*factor*Dy);  %Boundary Upper
   HyPMLC(1,:)=HyPMLC(2,:); %Boundary Left;

   %The following part is for ZONE D. 
   %H components.
   HxPMLD(:,2:NTy-1)=HxPMLD(:,2:NTy-1)-...
      Dt*(EzxPMLD(:,2:NTy-1)+EzyPMLD(:,2:NTy-1)-EzxPMLD(:,1:NTy-2)-EzyPMLD(:,1:NTy-2))/(mu0*Dy);
   
   HyPMLD(1:NPML-1,:)=exp(-Sigmax_yD(1:NPML-1,:)*Dt/e0).*HyPMLD(1:NPML-1,:)+...
      (1-exp(-Sigmax_yD(1:NPML-1,:)*Dt/e0))./(Sigmax_yD(1:NPML-1,:)*factor*Dx).*...
      (EzxPMLD(2:NPML,:)+EzyPMLD(2:NPML,:)-EzxPMLD(1:NPML-1,:)-EzyPMLD(1:NPML-1,:));
   
   HxPMLD(:,1)=expboundary*HxPMLD(:,1)-(1-expboundary)*...
      (EzxPMLD(:,1)+EzyPMLD(:,1)-EzxPML4(:,NPML)-EzyPML4(:,NPML))/(SigmaBound*factor*Dy);  %Boundary Down
   HxPMLD(:,NTy)=expboundary*HxPMLD(:,NTy)-(1-expboundary)*...
      (EzxPML2(:,1)+EzyPML2(:,1)-EzxPMLD(:,NTy-1)-EzyPMLD(:,NTy-1))/(SigmaBound*factor*Dy);  %Boundary Upper
   HyPMLD(NPML,:)=HyPMLD(NPML-1,:); %Boundary Right;

   %The following part is for ZONE 1.
   %H components.
   HxPML1(:,1:NPML-1)=exp(-Sigmay_x1(:,1:NPML-1)*Dt/e0).*HxPML1(:,1:NPML-1)-...
      (1-exp(-Sigmay_x1(:,1:NPML-1)*Dt/e0))./(Sigmay_x1(:,1:NPML-1)*factor*Dy).*...
      (EzxPML1(:,2:NPML)+EzyPML1(:,2:NPML)-EzxPML1(:,1:NPML-1)-EzyPML1(:,1:NPML-1));
   HxPML1(:,NPML)=HxPML1(:,NPML-1);
   
   HyPML1(2:NPML,:)=exp(-Sigmax_y1(2:NPML,:)*Dt/e0).*HyPML1(2:NPML,:)+...
      (1-exp(-Sigmax_y1(2:NPML,:)*Dt/e0))./(Sigmax_y1(2:NPML,:)*factor*Dx).*...
      (EzxPML1(2:NPML,:)+EzyPML1(2:NPML,:)-EzxPML1(1:NPML-1,:)-EzyPML1(1:NPML-1,:));
   HyPML1(1,:)=HyPML1(2,:);
   
   %The following part is for ZONE 2.
   %H components.
   HxPML2(:,1:NPML-1)=exp(-Sigmay_x2(:,1:NPML-1)*Dt/e0).*HxPML2(:,1:NPML-1)-...
      (1-exp(-Sigmay_x2(:,1:NPML-1)*Dt/e0))./(Sigmay_x2(:,1:NPML-1)*factor*Dy).*...
      (EzxPML2(:,2:NPML)+EzyPML2(:,2:NPML)-EzxPML2(:,1:NPML-1)-EzyPML2(:,1:NPML-1));
   HxPML2(:,NPML)=HxPML2(:,NPML-1);
   
   HyPML2(1:NPML-1,:)=exp(-Sigmax_y2(1:NPML-1,:)*Dt/e0).*HyPML2(1:NPML-1,:)+...
      (1-exp(-Sigmax_y2(1:NPML-1,:)*Dt/e0))./(Sigmax_y2(1:NPML-1,:)*factor*Dx).*...
      (EzxPML2(2:NPML,:)+EzyPML2(2:NPML,:)-EzxPML2(1:NPML-1,:)-EzyPML2(1:NPML-1,:));
   HyPML2(NPML,:)=HyPML2(NPML-1,:); %Boundary Right;
   
   %The following part is for ZONE 3.
   %H components.
   HyPML3(2:NPML,:)=exp(-Sigmax_y3(2:NPML,:)*Dt/e0).*HyPML3(2:NPML,:)+...
      (1-exp(-Sigmax_y3(2:NPML,:)*Dt/e0))./(Sigmax_y3(2:NPML,:)*factor*Dx).*...
      (EzxPML3(2:NPML,:)+EzyPML3(2:NPML,:)-EzxPML3(1:NPML-1,:)-EzyPML3(1:NPML-1,:));
   HyPML3(1,:)=HyPML3(2,:); %Boundary Left;

   HxPML3(:,2:NPML)=exp(-Sigmay_x3(:,2:NPML)*Dt/e0).*HxPML3(:,2:NPML)-...
      (1-exp(-Sigmay_x3(:,2:NPML)*Dt/e0))./(Sigmay_x3(:,2:NPML)*factor*Dy).*...
      (EzxPML3(:,2:NPML)+EzyPML3(:,2:NPML)-EzxPML3(:,1:NPML-1)-EzyPML3(:,1:NPML-1));
   HxPML3(:,1)=HxPML3(:,2); %Boundary Bottom;
   
   %The following part is for ZONE 4.
   %H components.
   HxPML4(:,2:NPML)=exp(-Sigmay_x4(:,2:NPML)*Dt/e0).*HxPML4(:,2:NPML)-...
      (1-exp(-Sigmay_x4(:,2:NPML)*Dt/e0))./(Sigmay_x4(:,2:NPML)*factor*Dy).*...
      (EzxPML4(:,2:NPML)+EzyPML4(:,2:NPML)-EzxPML4(:,1:NPML-1)-EzyPML4(:,1:NPML-1));
   HxPML4(:,1)=HxPML4(:,2); %Boundary Bottom;

   HyPML4(1:NPML-1,:)=exp(-Sigmax_y4(1:NPML-1,:)*Dt/e0).*HyPML4(1:NPML-1,:)+...
      (1-exp(-Sigmax_y4(1:NPML-1,:)*Dt/e0))./(Sigmax_y4(1:NPML-1,:)*factor*Dx).*...
      (EzxPML4(2:NPML,:)+EzyPML4(2:NPML,:)-EzxPML4(1:NPML-1,:)-EzyPML4(1:NPML-1,:));
   HyPML4(NPML,:)=HyPML4(NPML-1,:); %Boundary Right;

   %The following part is for inside region. 
   %Ez component.
   Ez=Ez+Dt*((Hy(2:NTx,1:NTy-1)-Hy(1:NTx-1,1:NTy-1))/Dx-...
             (Hx(1:NTx-1,2:NTy)-Hx(1:NTx-1,1:NTy-1))/Dy)./Ep; %Inside region.
          
   %The following part is for source
   Ez(1,Npy-(NMlat-1)/2:Npy+(NMlat-1)/2)=Ez(1,Npy-(NMlat-1)/2:Npy+(NMlat-1)/2)+sin(W*m*Dt);%*exp(-(m*W*Dt-3)^2); 
   %End of Source


   %The following part is for ZONE A. 
   %Ez component.
   EzxPMLA=EzxPMLA+Dt*(HyPMLA(2:NTx,:)-HyPMLA(1:NTx-1,:))/(e0*Dx);
   EzyPMLA(:,2:NPML)=exp(-Sigmay_zA(:,2:NPML)*Dt/e0).*EzyPMLA(:,2:NPML)-...
      (1-exp(-Sigmay_zA(:,2:NPML)*Dt/e0))./(Sigmay_zA(:,2:NPML)*Dy).*...
      (HxPMLA(:,2:NPML)-HxPMLA(:,1:NPML-1));    
   EzyPMLA(:,1)=exp(-Sigmay_zA(:,1)*Dt/e0).*EzyPMLA(:,1)-...
      (1-exp(-Sigmay_zA(:,1)*Dt/e0))./(Sigmay_zA(:,1)*Dy).*...
      (HxPMLA(:,1)-Hx(:,NTy));    
   
   %The following part is for ZONE B. 
   %Ez component.
   EzxPMLB=EzxPMLB+Dt*(HyPMLB(2:NTx,:)-HyPMLB(1:NTx-1,:))/(e0*Dx);
   EzyPMLB(:,1:NPML-1)=exp(-Sigmay_zB(:,1:NPML-1)*Dt/e0).*EzyPMLB(:,1:NPML-1)-...
      (1-exp(-Sigmay_zB(:,1:NPML-1)*Dt/e0))./(Sigmay_zB(:,1:NPML-1)*Dy).*...
      (HxPMLB(:,2:NPML)-HxPMLB(:,1:NPML-1));    
   EzyPMLB(:,NPML)=exp(-Sigmay_zB(:,NPML)*Dt/e0).*EzyPMLB(:,NPML)-...
      (1-exp(-Sigmay_zB(:,NPML)*Dt/e0))./(Sigmay_zB(:,NPML)*Dy).*...
      (Hx(:,1)-HxPMLB(:,NPML));    
   
   %The following part is for ZONE C. 
   %Ez component.
   EzxPMLC(1:NPML-1,:)=exp(-Sigmax_zC(1:NPML-1,:)*Dt/e0).*EzxPMLC(1:NPML-1,:)+...
      (1-exp(-Sigmax_zC(1:NPML-1,:)*Dt/e0))./(Sigmax_zC(1:NPML-1,:)*Dx).*...
      (HyPMLC(2:NPML,:)-HyPMLC(1:NPML-1,:));    
   EzxPMLC(NPML,:)=exp(-Sigmax_zC(NPML,:)*Dt/e0).*EzxPMLC(NPML,:)+...
      (1-exp(-Sigmax_zC(NPML,:)*Dt/e0))./(Sigmax_zC(NPML,:)*Dx).*...
      (Hy(1,:)-HyPMLC(NPML,:));    
   EzyPMLC=EzyPMLC-Dt*(HxPMLC(:,2:NTy)-HxPMLC(:,1:NTy-1))/(e0*Dy);

   %The following part is for ZONE D. 
   %Ez component.
   EzxPMLD(2:NPML,:)=exp(-Sigmax_zD(2:NPML,:)*Dt/e0).*EzxPMLD(2:NPML,:)+...
      (1-exp(-Sigmax_zD(2:NPML,:)*Dt/e0))./(Sigmax_zD(2:NPML,:)*Dx).*...
      (HyPMLD(2:NPML,:)-HyPMLD(1:NPML-1,:));    
   EzxPMLD(1,:)=exp(-Sigmax_zD(1,:)*Dt/e0).*EzxPMLD(1,:)+...
      (1-exp(-Sigmax_zD(1,:)*Dt/e0))./(Sigmax_zD(1,:)*Dx).*...
      (HyPMLD(1,:)-Hy(NTx,:));    
   EzyPMLD=EzyPMLD-Dt*(HxPMLD(:,2:NTy)-HxPMLD(:,1:NTy-1))/(e0*Dy);
   
   %The following part is for ZONE 1. 
   %Ez component.
   EzxPML1(1:NPML-1,:)=exp(-Sigmax_z1(1:NPML-1,:)*Dt/e0).*EzxPML1(1:NPML-1,:)+...
      (1-exp(-Sigmax_z1(1:NPML-1,:)*Dt/e0))./(Sigmax_z1(1:NPML-1,:)*Dx).*...
      (HyPML1(2:NPML,:)-HyPML1(1:NPML-1,:));    
   EzxPML1(NPML,:)=exp(-Sigmax_z1(NPML,:)*Dt/e0).*EzxPML1(NPML,:)+...
      (1-exp(-Sigmax_z1(NPML,:)*Dt/e0))./(Sigmax_z1(NPML,:)*Dx).*...
      (HyPMLA(1,:)-HyPML1(NPML,:)); 
   EzyPML1(:,2:NPML)=exp(-Sigmay_z1(:,2:NPML)*Dt/e0).*EzyPML1(:,2:NPML)-...
      (1-exp(-Sigmay_z1(:,2:NPML)*Dt/e0))./(Sigmay_z1(:,2:NPML)*Dy).*...
      (HxPML1(:,2:NPML)-HxPML1(:,1:NPML-1));    
   EzyPML1(:,1)=exp(-Sigmay_z1(:,1)*Dt/e0).*EzyPML1(:,1)-...
      (1-exp(-Sigmay_z1(:,1)*Dt/e0))./(Sigmay_z1(:,1)*Dy).*...
      (HxPML1(:,1)-HxPMLC(:,NTy));    
   
   %The following part is for ZONE 2. 
   %Ez component.
   EzxPML2(2:NPML,:)=exp(-Sigmax_z2(2:NPML,:)*Dt/e0).*EzxPML2(2:NPML,:)+...
      (1-exp(-Sigmax_z2(2:NPML,:)*Dt/e0))./(Sigmax_z2(2:NPML,:)*Dx).*...
      (HyPML2(2:NPML,:)-HyPML2(1:NPML-1,:));    
   EzxPML2(1,:)=exp(-Sigmax_z2(1,:)*Dt/e0).*EzxPML2(1,:)+...
      (1-exp(-Sigmax_z2(1,:)*Dt/e0))./(Sigmax_z2(1,:)*Dx).*...
      (HyPML2(1,:)-HyPMLA(NTx,:));    
   EzyPML2(:,2:NPML)=exp(-Sigmay_z2(:,2:NPML)*Dt/e0).*EzyPML2(:,2:NPML)-...
      (1-exp(-Sigmay_z2(:,2:NPML)*Dt/e0))./(Sigmay_z2(:,2:NPML)*Dy).*...
      (HxPML2(:,2:NPML)-HxPML2(:,1:NPML-1));    
   EzyPML2(:,1)=exp(-Sigmay_z2(:,1)*Dt/e0).*EzyPML2(:,1)-...
      (1-exp(-Sigmay_z2(:,1)*Dt/e0))./(Sigmay_z2(:,1)*Dy).*...
      (HxPML2(:,1)-HxPMLD(:,NTy));    
   
   %The following part is for ZONE 3.
   %Ez component.
   EzxPML3(1:NPML-1,:)=exp(-Sigmax_z3(1:NPML-1,:)*Dt/e0).*EzxPML3(1:NPML-1,:)+...
      (1-exp(-Sigmax_z3(1:NPML-1,:)*Dt/e0))./(Sigmax_z3(1:NPML-1,:)*Dx).*...
      (HyPML3(2:NPML,:)-HyPML3(1:NPML-1,:));    
   EzxPML3(NPML,:)=exp(-Sigmax_z3(NPML,:)*Dt/e0).*EzxPML3(NPML,:)+...
      (1-exp(-Sigmax_z3(NPML,:)*Dt/e0))./(Sigmax_z3(NPML,:)*Dx).*...
      (HyPMLB(1,:)-HyPML3(NPML,:));    

   EzyPML3(:,1:NPML-1)=exp(-Sigmay_z3(:,1:NPML-1)*Dt/e0).*EzyPML3(:,1:NPML-1)-...
      (1-exp(-Sigmay_z3(:,1:NPML-1)*Dt/e0))./(Sigmay_z3(:,1:NPML-1)*Dy).*...
      (HxPML3(:,2:NPML)-HxPML3(:,1:NPML-1));    
   EzyPML3(:,NPML)=exp(-Sigmay_z3(:,NPML)*Dt/e0).*EzyPML3(:,NPML)-...
      (1-exp(-Sigmay_z3(:,NPML)*Dt/e0))./(Sigmay_z3(:,NPML)*Dy).*...
      (HxPMLC(:,1)-HxPML3(:,NPML));    

   %The following part is for ZONE 4. 
   %Ez component.
   EzxPML4(2:NPML,:)=exp(-Sigmax_z4(2:NPML,:)*Dt/e0).*EzxPML4(2:NPML,:)+...
      (1-exp(-Sigmax_z4(2:NPML,:)*Dt/e0))./(Sigmax_z4(2:NPML,:)*Dx).*...
      (HyPML4(2:NPML,:)-HyPML4(1:NPML-1,:));    
   EzxPML4(1,:)=exp(-Sigmax_z4(1,:)*Dt/e0).*EzxPML4(1,:)+...
      (1-exp(-Sigmax_z4(1,:)*Dt/e0))./(Sigmax_z4(1,:)*Dx).*...
      (HyPML4(1,:)-HyPMLB(NTx,:));    
   EzyPML4(:,1:NPML-1)=exp(-Sigmay_z4(:,1:NPML-1)*Dt/e0).*EzyPML4(:,1:NPML-1)-...
      (1-exp(-Sigmay_z4(:,1:NPML-1)*Dt/e0))./(Sigmay_z4(:,1:NPML-1)*Dy).*...
      (HxPML4(:,2:NPML)-HxPML4(:,1:NPML-1));    
   EzyPML4(:,NPML)=exp(-Sigmay_z4(:,NPML)*Dt/e0).*EzyPML4(:,NPML)-...
      (1-exp(-Sigmay_z4(:,NPML)*Dt/e0))./(Sigmay_z4(:,NPML)*Dy).*...
      (HxPMLD(:,1)-HxPML4(:,NPML));    
   
   
   if mod(m,Meach)==0
      m
      if IsFigure==1
         %The following part is for visual figures. 
         EzAll=[EzxPML3+EzyPML3,EzxPMLC+EzyPMLC,EzxPML1+EzyPML1;
                EzxPMLB+EzyPMLB,Ez,EzxPMLA+EzyPMLA;
                EzxPML4+EzyPML4,EzxPMLD+EzyPMLD,EzxPML2+EzyPML2];
         HxAll=[HxPML3,HxPMLC,HxPML1;HxPMLB,Hx,HxPMLA;HxPML4,HxPMLD,HxPML2];
         HyAll=[HyPML3,HyPMLC,HyPML1;HyPMLB,Hy,HyPMLA;HyPML4,HyPMLD,HyPML2];

         figure(1);
         clf
         surf(X,Y,real(EzAll));
         caxis([-Colormax Colormax]);
         shading interp;
         axis([-0.5*a*MLatx-Dx*NPML 0.5*a*MLatx+Dx*NPML -0.5*a*MLaty-Dx*NPML 0.5*a*MLaty+Dx*NPML])
         zlabel('Ez');
         if IsMovie==1
            Mnum=Mnum+1;
            Movie(:,Mnum)=getframe;
         end
         view(0,90);
         axis off;
         pause(0.2)
      end
   end
   %toc
end

toc