How can I revise my 1D FDTD code

This is my code to simulate electromagnetic wave transmitting in a slab, and this is an example from CEM lecture at youtube, but code is not provided.

clc
close all
clear all

Nz = 94;
c0 = 299792458;
dz = 4.293e-3;
dt = 7.1599e-12;
L  = Nz*dz;
x  = linspace(0,L,Nz);

%Source parameters
STEPS = 4095;
f2    = 1.0e9;
NFREQ = 100;
FREQ  = linspace(0,f2,NFREQ);
tau   = 5e-10;
t0    = 3e-9;
nz_src = 3;

%Initial fourier transforms
K   = exp (-1i*2*pi*dt*FREQ);
REF = zeros(1,NFREQ);
TRN = zeros(1,NFREQ);
SRC = zeros(1,NFREQ);

%Initialize materials to free space
ER = ones(1, Nz);
UR = ones(1, Nz);
ER(13:83) = 6.0;
UR(13:83) = 2.0;

%Compute update coefficients
mEy = (c0*dt)./ER;
mHx = (c0*dt)./UR;

%Compute Gaussian source functions
t    = [0:STEPS-1]*dt;              %time axis
delt = 1.074e-11;                   %total delay between E and H
A    = -sqrt(ER(1)/UR(1));          %amplitude of H field
Esrc = exp(-((t-t0)/tau).^2);       %E field source
Hsrc = A*exp(-((t-t0+delt)/tau).^2);%H field source

% initialize fields
Ey=zeros(1,Nz);
Hx=zeros(1,Nz);

%Initialize boundary terms
H3=0; H2=0; H1=0;
E3=0; E2=0; E1=0;

%Main FDTD loop
for T = 1 : STEPS
    %Perform FDTD function
    % H update
    %Hx(1) = Hx(2);
    for i = 1:(Nz-1)
        Hx(i)=Hx(i)+mHx(i).*(Ey(i+1) - Ey(i))/dz;
    end
    Hx(Nz) = Hx(Nz) + mHx(Nz).*(E3 - Ey(Nz))/dz;
    Hx(nz_src-1)=Hx(nz_src-1)-mHx(nz_src-1).*Esrc(T)/dz; % TF/SF source
    H3=H2; H2=H1; H1=Hx(1); % ABC
    
    % E update
    %Ey(Nz) = Ey(Nz-1);
    for i = 2:Nz
        Ey(i) = Ey(i) + mEy(i).*(Hx(i) -Hx(i-1))/dz;
    end
    Ey(1) = Ey(1) + mEy(1)*(Hx(1) - H3)/dz;
    Ey(nz_src-1)=Ey(nz_src-1)-mEy(nz_src-1).*Hsrc(T)/dz; % TF/SF source
    E3=E2; E2=E1; E1=Ey(Nz); % ABC
    
    %Update Fourier Transforms 
    for nf = 1 : NFREQ 
     REF(nf) = REF(nf) + (K(nf)^T)*Ey(1); 
     TRN(nf) = TRN(nf) + (K(nf)^T)*Ey(Nz); 
     SRC(nf) = SRC(nf) + (K(nf)^T)*Esrc(T); 
    end 
    
    %Visualize
    figure(2)
    hold off
    plot(x,Ey,'b');
    axis([0 L -1 1]);
    grid on;
    hold on
    plot(x,Hx,'r');
    pause(0);
    T

end

%Finish fourier transforms
REF = REF*dt;
TRN = TRN*dt;
SRC = SRC*dt;

%Compute reflectance and transmittance
REF = abs(REF./SRC).^2;
TRN = abs(TRN./SRC).^2;
CON = REF + TRN;

figure(1)
plot(FREQ,REF)
hold on
plot(FREQ,TRN,'r-')
plot(FREQ,TRN+REF,'k-')
xlabel('Frequency(Hz)')
ylabel('Amplitude')
title('Impulse response of Etalon')
legend on
legend('Rx','Tx','Tx+Rx')
hold off
fh = figure(1);
set(fh, 'color', 'white');

There are some promblems with the boundary condition that I cannot fix. And when I change the value of nz_src, the simulation result will be totally different, the calue of E and H will reach to about 40! How could this even happen?

I referenced some code from matlab central.

Answer 1

[Sandeep Yadav Golla]

clc;
close all;
clear all;

Nz = 94;
c0 = 299792458;
dz = 4.293e-3;
dt = 7.1599e-12;
L  = Nz*dz;
x  = linspace(0,L,Nz);
%Source parameters
STEPS = 4095;
f2    = 1.0e9;
NFREQ = 100;
FREQ  = linspace(0,f2,NFREQ);
tau   = 5e-10;
t0    = 3e-9;
nz_src = 2;                            %%% corrected here %%%
%Initial fourier transforms
K   = exp (-1i*2*pi*dt*FREQ);
REF = zeros(1,NFREQ);
TRN = zeros(1,NFREQ);
SRC = zeros(1,NFREQ);
%Initialize materials to free space
ER = ones(1, Nz);
UR = ones(1, Nz);
ER(13:83) = 6.0;
UR(13:83) = 2.0;
%Compute update coefficients
mEy = (c0*dt)./ER;
mHx = (c0*dt)./UR;
%Compute Gaussian source functions
t    = [0:STEPS-1]*dt;              %time axis
delt = 1.074e-11;                   %total delay between E and H
A    = -sqrt(ER(1)/UR(1));          %amplitude of H field
Esrc = exp(-((t-t0)/tau).^2);       %E field source
Hsrc = A*exp(-((t-t0+delt)/tau).^2);%H field source
% initialize fields
Ey=zeros(1,Nz);
Hx=zeros(1,Nz);
%Initialize boundary terms
H2=0; H1=0;                        %%%% corrected here %%%
E2=0; E1=0;
%Main FDTD loop
for T = 1 : STEPS
    %Perform FDTD function
    % H update
    %Hx(1) = Hx(2);
    H2=H1; H1=Hx(1); % ABC
    for i = 1:(Nz-1)
        Hx(i)=Hx(i)+mHx(i)*(Ey(i+1) - Ey(i))/dz;
    end
    Hx(Nz) = Hx(Nz) + mHx(Nz)*(E2 - Ey(Nz))/dz;
    Hx(nz_src-1)=Hx(nz_src-1)-mHx(nz_src-1)*Esrc(T)/dz; % TF/SF source
    
    
    % E update
    %Ey(Nz) = Ey(Nz-1);
    E2=E1; E1=Ey(Nz); % ABC
    Ey(1) = Ey(1) + mEy(1)*(Hx(1) - H2)/dz;
    for i = 2:Nz
        Ey(i) = Ey(i) + mEy(i)*(Hx(i) -Hx(i-1))/dz;
    end
    Ey(nz_src)=Ey(nz_src)-mEy(nz_src)*Hsrc(T)/dz; % TF/SF source   %%% corrected here %%%
    
    
    %Update Fourier Transforms 
    for nf = 1 : NFREQ 
     REF(nf) = REF(nf) + (K(nf)^T)*Ey(1); 
     TRN(nf) = TRN(nf) + (K(nf)^T)*Ey(Nz); 
     SRC(nf) = SRC(nf) + (K(nf)^T)*Esrc(T); 
    end 
    
    %Visualize                       %%% wait some time  to see the simulations
    figure(2)
    hold off
    plot(x,Ey,'b');
    axis([0 L -1 1]);
    grid on;
    hold on
    plot(x,Hx,'r');
    pause(0);
end
%Finish fourier transforms
 REF = REF*dt;
 TRN = TRN*dt;
 SRC = SRC*dt;
%Compute reflectance and transmittance
REF = abs(REF./SRC).^2;
TRN = abs(TRN./SRC).^2;
CON = REF + TRN;
figure(1)
plot(FREQ,REF)
hold on
plot(FREQ,TRN,'r-')
plot(FREQ,TRN+REF,'k-')
xlabel('Frequency(Hz)')
ylabel('Amplitude')
title('Impulse response of Etalon')
legend on
legend('Rx','Tx','Tx+Rx')
hold off
fh = figure(1);
set(fh, 'color', 'white');

Answer 2

[alican kara]

https://www.youtube.com/watch?v=iugj1VC3QoU