function [V2] = ComputeIR(Dt, Dz, Nt, k2)
% [V2] = ComputeIR(Dt, Dz, Nt, k2)
  
  c0dt2 = (299792458 * Dt)^2;
  
  B = c0dt2 / (Dz * Dz);
  A = 2 - 2 * B - c0dt2 * k2;
  
  V2 = zeros(Nt+1,1);
  V3 = zeros(Nt+1,1);
  V4 = zeros(Nt+1,1);
  
  V2(1) = 0;
  V2(2) = B;
  V2(3) = A * B;
  V3(3) = B * B;
  
  for i = 4:Nt
    V2(i) = A * V2(i - 1) - V2(i - 2) + B * V3(i - 1);
    V3(i) = A * V3(i - 1) - V3(i - 2) + B * (V4(i - 1) + V2(i - 1));
    V4(i) = V2(i-2:-1:2)' * V3(3:i-1);
  end

  V2 = V2(2:end);
