Implementation of a TFSF source in four sides of FDTD

Hi all,

I am trying to correct my code by the implementation of the TFSF source in four sides of the FDTD problem, however, I have an issue in auxiliary griding of the code. I attached my code here, and would highly appreciate any help.

clear all; clc; close all;
% Grid Dimension in x (xdim) and y (ydim) directions
xN=400;
yN=xN;
%Total no of time steps
tN=600;
%Position of the source (center of the domain)
x_s=100;
y_s=100;
%Courant stability factor
S=1/(2^0.5);
% S=1/sqrt(3);
% Parameters of free space (permittivity and permeability and speed of
% light) are all not 1 and are given real values
epsilon0=(1/(36*pi))*1e-9;
mu0=4*pi*1e-7;
c=3e+8;
frequency=5e+13;
lambda=c/frequency;
% Spatial grid step length (spatial grid step= 1 micron and can be changed)
delx=lambda/10;
% Temporal grid step obtained using Courant condition
delt=S*delx/c;
% Initialization of field matrices
Ez=zeros(xN,yN);
Hy=zeros(xN,yN);
Hx=zeros(xN,yN);

Ezf=zeros(xN,yN);
Hyf=zeros(xN,yN);
Hxf=zeros(xN,yN);

Ezinc=zeros(xN,yN);
Hyinc=zeros(xN,yN);
Hxinc=zeros(xN,yN);

% Initialization of permittivity and permeability matrices
epsilon=epsilon0*ones(xN,yN);
mu=mu0*ones(xN,yN);
% Initializing electric and magnetic conductivity matrices
sigma=4e-4*ones(xN,yN);
sigma_star=4e-4*ones(xN,yN);
%Choice of nature of source
gaussian=0;
sine=0;
% The user can give a frequency of his choice for sinusoidal (if sine=1 above) waves in Hz 

impulse=0;
%Choose any one as 1 and rest as 0. Default (when all are 0): Unit time step
%Multiplication factor matrices for H matrix update to avoid being calculated many times 
%in the time update loop so as to increase computation speed
A=((mu-0.5*delt*sigma_star)./(mu+0.5*delt*sigma_star)); 
B=(delt/delx)./(mu+0.5*delt*sigma_star);
                          
%Multiplication factor matrices for E matrix update to avoid being calculated many times 
%in the time update loop so as to increase computation speed                          
C=((epsilon-0.5*delt*sigma)./(epsilon+0.5*delt*sigma)); 
D=(delt/delx)./(epsilon+0.5*delt*sigma);
F=(delt./epsilon)./(1+0.5*sigma*delt./epsilon)+ones(xN,yN); % Jsource coeffecient
Cb=(c*delt-delx)/(c*delt+delx);
one_way=1; %1 for One way BC, 0 for PEC
jsource=0; %1=Jsource ON, 0=Jsource OFF

%!-----------------Source parameters---------------------------------%
t=(1:tN)*delt;
Zs=50; % position index of source
delta=20*delt;

% Update loop begins
for n=1:1:tN

    if jsource==1
        % J source input
        if n<=30
            f(x_s,y_s)=sin(2*pi*frequency*(n)*delt)*exp(-(n-30)^2/30^2);
        else
            f(x_s,y_s)=sin(2*pi*frequency*(n)*delt);
        end
        end
    %if source is impulse or unit-time step 
    if impulse==1 
        if n==1
        Ez(x_s,y_s)=1;
        else
        Ez(x_s,y_s)=0;    
        end
    end
    
    %% Update fields
    %Vector update instead of for-loop for Hy and Hx fields
    
    Ez(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
    Hy(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hy(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ez(1+1:xN-1+1,1:yN-1)-Ez(1:xN-1,1:yN-1));
    Hx(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hx(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ez(1:xN-1,1+1:yN-1+1)-Ez(1:xN-1,1:yN-1));
    Ez(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ez(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hy(1+1:xN-1+1,1+1:yN-1+1)-Hy(1:xN-1,1+1:yN-1+1)-Hx(1+1:xN-1+1,1+1:yN-1+1)+Hx(1+1:xN-1+1,1:yN-1));
    
    
%      Ezinc(Zs,xN/2-50:xN/2+50)=cos(frequency*delt*(n-3));
     Ezinc(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
%      Ezinc(Zs+1,Zs+2:yN-Zs-2)=(1-exp(-((delt*n).^2)/(delta^2)))*cos(frequency*delt*(n-3));
%    Hyinc(Zs,:)=Am*cos(frequency.*(delt*n+s));
     Hyinc(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hyinc(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ezinc(1+1:xN-1+1,1:yN-1)-Ezinc(1:xN-1,1:yN-1));
     Hxinc(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hxinc(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ezinc(1:xN-1,1+1:yN-1+1)-Ezinc(1:xN-1,1:yN-1));
     Ezinc(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ezinc(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hyinc(1+1:xN-1+1,1+1:yN-1+1)-Hyinc(1:xN-1,1+1:yN-1+1)-Hxinc(1+1:xN-1+1,1+1:yN-1+1)+Hxinc(1+1:xN-1+1,1:yN-1));
    
    
    Hy(Zs-1,Zs:xN-Zs)=Hy(Zs-1,Zs:xN-Zs)-B(Zs,Zs:xN-Zs).*Ezinc(Zs,Zs:xN-Zs);  
    Hy(xN-Zs+1,Zs:xN-Zs)=Hy(xN-Zs+1,Zs:xN-Zs)+B(xN-Zs,Zs:xN-Zs).*Ezinc(xN-Zs,Zs:xN-Zs);  
    Hx(Zs:xN-Zs,Zs-1)=Hx(Zs:xN-Zs,Zs-1)+B(Zs:xN-Zs,Zs).*Ezinc(Zs:xN-Zs,Zs);  
    Hx(Zs:xN-Zs,xN-Zs+1)=Hx(Zs:xN-Zs,xN-Zs+1)-B(Zs:xN-Zs,xN-Zs).*Ezinc(Zs:xN-Zs,xN-Zs); 
      
    Ez(Zs,Zs:xN-Zs)=Ez(Zs,Zs:xN-Zs)-D(Zs,Zs:xN-Zs).*Hyinc(Zs-1,Zs:xN-Zs);
    Ez(xN-Zs,Zs:xN-Zs)=Ez(xN-Zs,Zs:xN-Zs)+D(xN-Zs-1,Zs:xN-Zs).*Hyinc(xN-Zs-1,Zs:xN-Zs);
    Ez(Zs:xN-Zs,Zs)=Ez(Zs:xN-Zs,Zs)+D(Zs:xN-Zs,Zs).*Hxinc(Zs:xN-Zs,Zs-1);
    Ez(Zs:xN-Zs,xN-Zs)=Ez(Zs:xN-Zs,xN-Zs)-D(Zs:xN-Zs,xN-Zs).*Hxinc(Zs:xN-Zs,xN-Zs+1);
     
    if jsource ==1
        Ez(x_s,y_s)=Ez(x_s,y_s)-F(x_s,y_s)*f(x_s,y_s); %J_source update
    end
    
    %% temporary storage of fields for getting Ez(n+1) (required for One way BC)
    Hyf=Hy;
    Hxf=Hx;
    Ezf=Ez;
    Hyf(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hyf(1:xN-1,1:yN-1)+B(1:xN-1,1:yN-1).*(Ezf(1+1:xN-1+1,1:yN-1)-Ezf(1:xN-1,1:yN-1));
    Hxf(1:xN-1,1:yN-1)=A(1:xN-1,1:yN-1).*Hxf(1:xN-1,1:yN-1)-B(1:xN-1,1:yN-1).*(Ezf(1:xN-1,1+1:yN-1+1)-Ezf(1:xN-1,1:yN-1));
    Ezf(1+1:xN-1+1,1+1:yN-1+1)=C(1+1:xN-1+1,1+1:yN-1+1).*Ezf(1+1:xN-1+1,1+1:yN-1+1)+D(1+1:xN-1+1,1+1:yN-1+1).*(Hyf(1+1:xN-1+1,1+1:yN-1+1)-Hyf(1:xN-1,1+1:yN-1+1)-Hxf(1+1:xN-1+1,1+1:yN-1+1)+Hxf(1+1:xN-1+1,1:yN-1));
    
    if jsource==1
        Ezf(x_s,y_s)=Ezf(x_s,y_s)-F(x_s,y_s)*f(x_s,y_s);
    end

    %% Boundary condition
if one_way==1
    Ez(1,:)=Ez(2,:)+Cb.*(Ezf(2,:)-Ez(1,:));             %at x=0
    Ez(xN,:)=Ez(xN-1,:)+Cb.*(Ezf(xN-1,:)-Ez(xN,:));     %at x=l
    Ez(:,1)=Ez(:,2)+Cb.*(Ezf(:,2)-Ez(:,1));             %at y=0
    Ez(:,yN)=Ez(:,yN-1)+Cb.*(Ezf(:,yN-1)-Ez(:,yN));     %at y=l

else
    % Perfect Electric Conductor boundary condition
       Ez(1,:)=0;
       Ez(xN,:)=0;
       Ez(:,1)=0;
       Ez(:,yN)=0;  
end
    %% Hard Source type
    if impulse==0 && jsource==0
        % If unit-time step
%         if gaussian==0 && sine==0
%             Ez(x_s,y_s)=0;
%         end
        %if sine
        if sine==1
            tstart=1;
            N_lambda=c/(frequency*delx);
             Ez(x_s,y_s)=sin(((2*pi*(c/(delx*N_lambda))*(n-tstart)*delt)));
%              Ez(:,y_s)=sin(((2*pi*(c/(delx*N_lambda))*(n-tstart)*delt)));
        end
        %if gaussian
        if gaussian==1
            if n<=42
                Ez(x_s,y_s)=(10-15*cos(n*pi/20)+6*cos(2*n*pi/20)-cos(3*n*pi/20))/32;
            else
                Ez(x_s,y_s)=0;
            end
        end
        
    end
    
    
    %% Movie plot of Ez
    if jsource==1
    imagesc(delx*(1:1:xN)*1e+6,(1e+6*delx*(1:1:yN))',Ez',[-0.05 0.05]); colormap jet;
    else
    imagesc(delx*(1:1:xN)*1e+6,(1e+6*delx*(1:1:yN))',Ez',[-1 1]); colormap jet;colorbar
    end

    title(['\fontsize{20}E_{z} at t = ',num2str(round(n*delt*1e+15)),' fs']); 
    xlabel('x (in um)','FontSize',20);
    ylabel('y (in um)','FontSize',20);
    set(gca,'FontSize',20);
    getframe;
end