Successive FFT2 transforms cause an increase in amplitude even though each individual transform is normalised.

I will be very grateful if anyone can help here. I have some code to simulate Fourier Optics and need to achieve this with a series of Fourier transforms over a short distance where I will eventually take account of turbulence at each step.

Unfortunately I have noticed that when I perform repeated FFT2 transforms on my gaussian beam profile, the resulting beam's amplitude if greatly increased even tho I am normalising the transformed function each step of the way.

My code is as follows:

 %%Fourier Optics Code: ===================================================

 clear;                  % clear workspace
 steps = 2^8;            % Resolution (256)
 lim = 1;                % Field of View (m)
 scale = 2/steps;        % Normalisation factor for fft & ifft
 ii = sqrt(-1);          % Neat Imaginary unit

 % Set Up Co-ordinate Grids: -----------------------------------------------
 % Direct Space Window 2D:
 x = linspace(-lim, lim, steps); y = x; [X,Y]=meshgrid(x,y);
 Xs = mean(diff(x));     % Sampling Interval  (delta x)
 Ks = 1/Xs;              % Sampling Frequency (per distance)
 Kn = Ks/2;              % Nyquist Frequency 

 % K-Space Window 2D: 
 kx = linspace(-1, 1, steps)*Kn; ky = kx; [KX,KY] = meshgrid(kx,ky);

 % Beam Properties: --------------------------------------------------------
 lambda = 900E-06;       % Wavelenght (meters)
 w0 = 0.01;              % Beam Wasit (meters)
 k0 = 2*pi/lambda;       % Wavenumber (meters^-1)
 z0 = (k0*w0^2)/2;       % Rayleigh Range

 % Propogation Distance: ---------------------------------------------------
 dz = 5;                 % meters
 zsteps = 100;           % number of propagation steps taken

 % Neat Veriable for Cylindrical Polars: -----------------------------------
 R = sqrt(X.^2+Y.^2);    % Radial Magnitude
 K = sqrt(KX.^2+KY.^2);  % "Radial" K-vector
 dGouy = -atan(dz/z0);   % Gouy Phase progression due to dz propagation

 % Making Beam: ------------------------------------------------------------
 func = exp((-1)*(R./w0).^2);    % Gaissian Beam (Waist = w0 in plane z = 0)

 % Propogate Beam ---------------------------------------------------------
 next = func; peaks = zeros(length(zsteps),1);
 for k = 1:zsteps   
      if k == 1;  % Ensures Beam in centre of Field of View
           Func = fftshift(fft2(next)*scale^2);
      else
           Func = fft2(next)*scale^2;
      end

      % Fourier Optics:
      H = exp(ii*k0*dz) .* exp(-ii*(dz.*K.^2)./(2*k0));    % Propagator 
      HF = H .* Func;                                 
      next = (ifft2(HF)/scale^2)*dGouy;  
      peaks(k) = (abs(max(max(next))));   % Record the Max of the output
end

% Plotting =================================================================
fig = figure(1);

subplot(1,3,1); 
surf(X,Y,abs(func),'EdgeColor','none');
title('Original Beam');
xlabel('Direct Space FOV (m)');
ylabel('Direct Space FOV (m)');
axis('tight');
view([0,0,1]);

subplot(1,3,2);
surf(X,Y,abs(next),'EdgeColor','none');
title('Propagated Beam');
xlabel('Direct Space FOV (m)');
ylabel('Direct Space FOV (m)');
axis('tight');

subplot(1,3,3);
plot(diff(peaks))
title('Increase in Max Amplitude');
xlabel('Number of Steps in Z');
ylabel('Change in Amplitude');
axis('tight');

which outputs something like this:

As you can see the amplitude of the beam is increasing exponentially with propagation which isn't right! The amplitude should decrease with distance, and the increase shouldn't be from the scaling of the FFT2 and IFFT2 functions as I have tested my scaling method extensively in other programs. If any one is able to shed any light on this issue or notice something I may be doing wrong, I would be extraordinarily greatful!