Connective boundary of 2D FDTD PML

%100x100

iblcb=11;             % array index of connective boundary, left %
ibrcb=ib-10;                                       % right %
jbfcb=11;                                           % front %
jbbcb=jb-10;                                      % back %

for n=1:7.0*tau
 source(n)=10*(-(2/tau^2)*n+6/tau)*exp(-((n-delay)^2/tau^2));  %gauss%
end

ez(12,50)=source(n);

% connective boundary  with Hx & Hy

    hx(iblcb:ibrcb,jbfcb)=hx(iblcb:ibrcb,jbfcb)...             
       -dt/muz*((ez(iblcb:ibrcb,jbfcb)-ez(iblcb:ibrcb,jbfcb-1))/dx)+dt/muz/dx*ez(iblcb:ibrcb,jbfcb);                     % front CB %
    
    hx(iblcb:ibrcb,jbbcb+1)=hx(iblcb:ibrcb,jbbcb+1)...         
       -dt/muz*((ez(iblcb:ibrcb,jbbcb+1)-ez(iblcb:ibrcb,jbbcb))/dx)-dt/muz/dx*ez(iblcb:ibrcb,jbbcb);                    % back CB %
    
    hy(iblcb,jbfcb:jbbcb)=hy(iblcb,jbfcb:jbbcb)...              
       +dt/muz*((ez(iblcb,jbfcb:jbbcb)-ez(iblcb-1,jbfcb:jbbcb))/dx)-dt/muz/dx*ez(iblcb,jbfcb:jbbcb);                   % left CB %
    
    hy(ibrcb+1,jbfcb:jbbcb)=hy(ibrcb+1,jbfcb:jbbcb)...     
        +dt/muz*((ez(ibrcb+1,jbfcb:jbbcb)-ez(ibrcb,jbfcb:jbbcb))/dx)+dt/muz/dx*ez(ibrcb,jbfcb:jbbcb);                   % right CB %



% connective boundary   with  Ez

    ez(iblcb,jbfcb:jbbcb)=ez(iblcb,jbfcb:jbbcb)+dt/epsz*((hy(iblcb+1,jbfcb:jbbcb)-hy(iblcb,jbfcb:jbbcb))/dx...
             -(hx(iblcb,jbfcb+1:jbbcb+1)-hx(iblcb,jbfcb:jbbcb))/dx)-(dt/epsz/dx)*hy(iblcb,jbfcb:jbbcb);                                  % left CB %
           
    ez(ibrcb,jbfcb:jbbcb)=ez(ibrcb,jbfcb:jbbcb)+dt/epsz*((hy(ibrcb+1,jbfcb:jbbcb)-hy(ibrcb,jbfcb:jbbcb))/dx...
             -(hx(ibrcb,jbfcb+1:jbbcb+1)-hx(ibrcb,jbfcb:jbbcb))/dx)+(dt/epsz/dx)*hy(ibrcb+1,jbfcb:jbbcb);                             % right CB %
          
    ez(iblcb:ibrcb,jbfcb)=ez(iblcb:ibrcb,jbfcb)+dt/epsz*((hy(iblcb+1:ibrcb+1,jbfcb)-hy(iblcb:ibrcb,jbfcb))/dx...
             -(hx(iblcb:ibrcb,jbfcb+1)-hx(iblcb:ibrcb,jbfcb))/dx)+(dt/epsz/dx)*hx(iblcb:ibrcb,jbfcb);                                  % front CB %
          
    ez(iblcb:ibrcb,jbbcb)=ez(iblcb:ibrcb,jbbcb)+dt/epsz*((hy(iblcb+1:ibrcb+1,jbbcb)-hy(iblcb:ibrcb,jbbcb))/dx...
             -(hx(iblcb:ibrcb,jbbcb+1)-hx(iblcb:ibrcb,jbbcb))/dx)-(dt/epsz/dx)*hx(iblcb:ibrcb,jbbcb+1);                              % back CB %


% 4 corners

    ez(iblcb,jbfcb)=ez(iblcb,jbfcb)+dt/epsz*((hy(iblcb+1,jbfcb)-hy(iblcb,jbfcb))/dx...
       -(hx(iblcb,jbfcb+1)-hx(iblcb,jbfcb))/dx)-(dt/epsz/dx)*hy(iblcb,jbfcb)+(dt/epsz/dx)*hx(iblcb,jbfcb);
    
    ez(iblcb,jbbcb)=ez(iblcb,jbbcb)+dt/epsz*((hy(iblcb+1,jbbcb)-hy(iblcb,jbbcb))/dx...
       -(hx(iblcb,jbbcb+1)-hx(iblcb,jbbcb))/dx)-(dt/epsz/dx)*hy(iblcb,jbbcb)-(dt/epsz/dx)*hx(iblcb,jbbcb+1);
    
    ez(ibrcb,jbbcb)=ez(ibrcb,jbbcb)+dt/epsz*((hy(ibrcb+1,jbbcb)-hy(ibrcb,jbbcb))/dx...
       -(hx(ibrcb,jbbcb+1)-hx(ibrcb,jbbcb))/dx)+(dt/epsz/dx)*hy(ibrcb+1,jbbcb)-(dt/epsz/dx)*hx(ibrcb,jbbcb+1);
    
    ez(ibrcb,jbfcb)=ez(ibrcb,jbfcb)+dt/epsz*((hy(ibrcb+1,jbfcb)-hy(ibrcb,jbfcb))/dx...
       -(hx(ibrcb,jbfcb+1)-hx(ibrcb,jbfcb))/dx)+(dt/epsz/dx)*hy(ibrcb+1,jbfcb)+(dt/epsz/dx)*hx(ibrcb,jbfcb);