I am trying to simulate the propagation of a gaussian beam through a lens using an FFT approach.
I tried to implement the approach described by Couairon in this paper at page 43: https://link.springer.com/article/10.1140/epjst/e2011-01503-3
Here is a description of the algorithm: 
.
However, when I propagate the beam from a lens to the focal spot, the intensity I get is much less than if the theoretical intensity that I can compute using the theoretical formulas(https://en.wikipedia.org/wiki/Gaussian_beam):
$$ I(r,z)=I_0 (\frac{w_0}{w(z)})^2 e^{(-2r^2/w^2(z)}$$
where $I_0$ is the intensity at the center of the beam at its waist.
Here is the code I wrote:
Wavelength=532e-9
n0=1.0 #Freespace index of refraction
f=100e-3 # Lens focal length
r_max=10e-3 # Limits of the grid
w_lens=1e-3 #Waist size at the lens
N_r=2**12 # Number of points
r_grid = np.linspace(-r_max, r_max, N_r) # Spatial grid
k0 = (n0) * 2 * np.pi / Wavelength # Freespace wavevector
dr = abs(r_grid[1] - r_grid[0])
E_initial= 1* np.exp(-(r_grid ** 2) / (w_lens ** 2) - 1j * k0 * r_grid ** 2 / (2 * f))
dz=f # Distance of propagation (propagating to the focal spot)
delta = dz / (4 * k0 * dr ** 2)
# Steps described in last item of step 2
kjdr1 = [2 * np.pi * j / N_r for j in np.arange(0, N_r / 2)]
kjdr2 = [2 * np.pi * (-1 + j / N_r) for j in np.arange(N_r / 2, N_r)]
kjdr = np.concatenate([kjdr1, kjdr2])
A = np.exp(-2 * 1j * delta * (kjdr) ** 2)
E_F = np.fft.fft(np.fft.fftshift(E_initial),norm="ortho")
E_F = (E_F * A) # Convolution
E_Final = np.fft.ifftshift(np.fft.ifft(E_F,norm="ortho"))
If you use these parameters, the ratio $\frac{I_{focus}}{I_{lens}}$ is 56.64 . However, if you use the theoretical relation, you get 3488.20.
Do you see anything wrong in my implementation of the algorithm? Anything you see that could explain why the ratio $\frac{I_{focus}}{I_{lens}}$ seems to be way off?
PS:My code uses field E but I convert E to I using $I=2*n_0*\epsilon_0*c*|E|^2$
