Computing a Hankel transform pair

I'm attempting to reproduce figure 1b from the following paper http://preprints.acmac.uoc.gr/33/1/ by utilising the hankel transform pair they provided in their published paper:

I have been able to successfully code (2) and produce the spectrum plot as in figure 1d, however i'm having difficulty correctly reproducing the intensity profile 1b, i've been trying for a while now and am currently stuck :/ i'm not really sure what i've done wrong, i have a feeling it's the way my loop functions when substituting the values of u(k) into the integral but im unsure how i'd compute (1) any other way.

attached below is my code:

 clc;       % clears command window
 clear all; % clear all variables from memory
 close all; % close all figure windows

 %%------------------------------------------------------------------------
 % Forward Hankel transform finding u(k)

 r_0 = 10;               % initial radius of the main ring
 alpha = 0.05;           % decay constant
 N = 1500;               % sets number of points
 k1=linspace(0,8,N);     % sets arbitrary interval for k

 % computes integral
 I1 = zeros(size(k1));   % pre-allocate matrix for integral values
 u0 = @(r) airy(r_0-r).*exp(alpha.*(r_0 - r));       % Radially symmetric Airy      beam
 f1 = @(r,k) 2*pi .* r .* besselj(0,k.*r).* u0(r);   % Contents of integral

 for i = 1:length(k1);
    temp = @(r) f1(r,k1(i));
    I1(i) = quadgk(temp,0,Inf);
 end

 plot(k1,I1); %plots the specturm
 ylabel('Spectrum'); xlabel('k');

 %%------------------------------------------------------------------------
 % Inverse Hankel transform finding u(r,z)

 f2 = @(x,r,z,I1) (1./(2.*pi)) .* x .* besselj(0,x.*r).* exp((-1i .* x.^2) .*  z./2) .* I1;

 r = linspace(-20,20,N);
 z = linspace(0,10,N);

 I2 = zeros(size(r));
 for k = 1:length(r)
   temp = @(x) f2(x,r(k),z(k),I1(k)); 
   I2(k) = quadgk(temp,0,50); %limit is 50 as integral doesnt converge
 end

 % plots the intensity profile of the wave
 plot(z,abs(I2).^2); ylabel('Intensity');xlabel('z');
 end