Calculation of reflection coefficient using FFT in FDTD code (matlab)

[kotsoss]

I want to obtain the reflection coefficient using FFT in a transmission line, the source was excited with a Gaussian pulse. Τhis is an abstract from a paper that I have as guide:

"In order to obtain the reflection coefficient, the coaxial cable was excited with a Gaussian pulse of the form
Vs (t) = exp(-(t − 3Ts)/Ts)
Ts =1/2fmax
which was added at every time step to the voltage node located at z = 6.340 cm. The number of time steps as well as the step size
that were used are given in the first row of Table III. The voltage waveform was sampled at the node located at z = 8.622 cm
in order to compute the reflection coefficient Γnum from a discrete Fourier transform of the incident and reflected pulses."


I have calculated the reflection coefficient but I think the results are wrong. Anyone knows whats wrongs? Do you have any code that I can see? It would help me very much.
So here is my code:

st=100;                  %trasmission line of 100 length
i=zeros(1,st);
v=zeros(1,st);


dz=0.01;
l=250;
dt=dz*sqrt(l*c);

for k=1:49         
c(k)=0.1;
end

for k=50:st
c(k)=0.3;
end

%matrixes for the storaging the values of the incident and reflected waves.
vr=zeros(1,18);            
vi=zeros(1,18);
vin=zeros(1,18);
vref=zeros(1,18);
sun=zeros(1,18);
r=0;
y=0;

for n=1:97
   
   t=n*dt;

  
    for k=1:st-1                             
        v(k)=v(k)-dt/(c(k)*dz)*(i(k+1)-i(k));      %calculate the voltage
    end
    
     
   
    v(st)=0;          
   fs=exp( -  ((t-1)/(0.3))^2   );         %gaussian pulse  
   v(1)=fs;
   
    for k=2:st
        i(k)=i(k)-dt/(l*dz)*(v(k)-v(k-1));   %calculate the current
    end

    if n>=23 && n<=40   %the time period where the incedent wave pass from the node k=21 
      r=r+1;                       
      vi(r)=v(21);          %save the values of the voltage of the incedent wave in a matrix
      vin(r)=fft(vi(r));     %calculate the fft of the incedent wave 
   end

   if n>=80 && n<=97       %the time period where the reflected wave pass from the node k=21 
      y=y+1;
      vr(y)=v(21);              %save the values of the voltage of the reflected wave in a matrix
      vref(y)=fft(vr(y));       %calculate the fft of the reflected wave
    end
   
    %plot(v);
    axis([0 st -2 2]);
    title(['Time = ',num2str(n)]);
    frame = getframe;
   
end

sun=abs(vref./vin);         %calculate the reflection coefficient
plot(sun)