How to implement N number of realizations of a for loop?

The for loop in the attached code runs only once. I need the simulation to execute N times. How should I modify the code to achieve this?

I need the Uout as a 51*51 matrix.

% vChGB Sine Hyperbolic Gaussian Beam
clc;close all;clear all;
lambda = 532; % wavelength in nanometers
k = 2*pi/lambda; % wave vector
z=[1.20 	2.014	4.91	6.638	8.097	10.179	12.804	14.491	16.594	18.136	20.468	22.842	26.091	28.695	31.132	32.486	34.735	36.173	38.068	40.256	42.089	44.672	46.358	48.066	50.045	52.816	54.523	56.191	58.919	61.127	62.085	64.106	66.127	67.856	70.084	72.272	74.189	76.063	78.064	80.522	82.876	85.856	86.147	88.189	90.606	92.606	94.584	96.23	98.336	99.794  100									];
z_diff = (diff(z)) ; 
z_diff = [0 z_diff]; 
ti=1  

ti = 1

beta=(pi/180)*ti ;  
Lx = .1 ;Ly = Lx;             
Nx = 51  ;Ny = Nx ;          
dx = Lx/Nx;dy = Ly/Ny;          
x = (-Nx/2:Nx/2-1)*dx;
y = (-Ny/2:Ny/2-1)*dy;
kx =  (-Nx/2:Nx/2-1)/(Nx.*dx);
ky =  (-Ny/2:Ny/2-1)/(Ny.*dy);
kz = ones(1,Nx);
k0 = sqrt(kx.^2 + ky.^2 + kz.^2  );
Kx=cos(beta)*kx + sin(beta)*kz;
Ky=ky;
Kz=-sin(beta)*kx + cos(beta)*kz;
delfx = 1/(Lx);%frequncy shift
delfy = delfx;
Phi = 0.0001; 

F=2*(pi^2)*(k^2).*z_diff.*Phi;
x0 = linspace(-.5,.5,Nx);y0=x0;[X0,Y0] = meshgrid(x0,y0);
w0 = 0.05;
Uin=exp(-(X0.^2+Y0.^2)./w0.^2) ;
x = x + z_diff(1)*tan(beta);
R = randn(Nx,Ny) ;% Gaussian random phase      
R = R - mean(R(:,:));% zero mean
R = R/std(R(:,:));   % unit variance
screen1 = exp(1i*((Kx.*x) + (Ky.*y))).*delfx.*delfy.*sqrt(F(1)).*R ; 
Uout = Uin.*exp(1i.*screen1);
c=0.7;
Ri=z_diff./cos(beta);
screens = zeros(Nx, Ny, length(F)-1);   
I_array = zeros(Nx, Ny, length(F)-1);   
I_array(:,:,1) = Uout.*conj(Uout);
for i = 2:length(F)
        x = x + z_diff(i)*tan(beta);
        screens(:,:,i) = screens(:,:,i-1) + exp(1i*((kx.*x)+(Ky.*y))).*delfx.*delfy.*sqrt(F(i)).*R;
        Hi = exp(1i*Kz*Ri(i) - c.*Ri(i)); 
        Ui1 = fftshift(fft2(Uout));     
        Ui2 = Hi.*Ui1;                
        Ui3 = ifft2(ifftshift(Ui2));     
        Uio = Ui3.*exp(1i.*screens(:,:,i));  
        I = Uio.*conj(Uio);              
        I_array(:,:,i) = I;            
        Uout = Uio;  
end

[Matt J]

The output below shows that the loop already runs 50 times. Not sure why you think otherwise.

% vChGB Sine Hyperbolic Gaussian Beam
clc;close all;clear all;
lambda = 532; % wavelength in nanometers
k = 2*pi/lambda; % wave vector
z=[1.20 	2.014	4.91	6.638	8.097	10.179	12.804	14.491	16.594	18.136	20.468	22.842	26.091	28.695	31.132	32.486	34.735	36.173	38.068	40.256	42.089	44.672	46.358	48.066	50.045	52.816	54.523	56.191	58.919	61.127	62.085	64.106	66.127	67.856	70.084	72.272	74.189	76.063	78.064	80.522	82.876	85.856	86.147	88.189	90.606	92.606	94.584	96.23	98.336	99.794  100									];
z_diff = (diff(z)) ; 
z_diff = [0 z_diff]; 
ti=1  

ti = 1

beta=(pi/180)*ti ;  
Lx = .1 ;Ly = Lx;             
Nx = 51  ;Ny = Nx ;          
dx = Lx/Nx;dy = Ly/Ny;          
x = (-Nx/2:Nx/2-1)*dx;
y = (-Ny/2:Ny/2-1)*dy;
kx =  (-Nx/2:Nx/2-1)/(Nx.*dx);
ky =  (-Ny/2:Ny/2-1)/(Ny.*dy);
kz = ones(1,Nx);
k0 = sqrt(kx.^2 + ky.^2 + kz.^2  );
Kx=cos(beta)*kx + sin(beta)*kz;
Ky=ky;
Kz=-sin(beta)*kx + cos(beta)*kz;
delfx = 1/(Lx);%frequncy shift
delfy = delfx;
Phi = 0.0001; 

F=2*(pi^2)*(k^2).*z_diff.*Phi;
x0 = linspace(-.5,.5,Nx);y0=x0;[X0,Y0] = meshgrid(x0,y0);
w0 = 0.05;
Uin=exp(-(X0.^2+Y0.^2)./w0.^2) ;
x = x + z_diff(1)*tan(beta);
R = randn(Nx,Ny) ;% Gaussian random phase      
R = R - mean(R(:,:));% zero mean
R = R/std(R(:,:));   % unit variance
screen1 = exp(1i*((Kx.*x) + (Ky.*y))).*delfx.*delfy.*sqrt(F(1)).*R ; 
Uout = Uin.*exp(1i.*screen1);
c=0.7;
Ri=z_diff./cos(beta);
screens = zeros(Nx, Ny, length(F)-1);   
I_array = zeros(Nx, Ny, length(F)-1);   
I_array(:,:,1) = Uout.*conj(Uout);
for i = 2:length(F)
    i
        x = x + z_diff(i)*tan(beta);
        screens(:,:,i) = screens(:,:,i-1) + exp(1i*((kx.*x)+(Ky.*y))).*delfx.*delfy.*sqrt(F(i)).*R;
        Hi = exp(1i*Kz*Ri(i) - c.*Ri(i)); 
        Ui1 = fftshift(fft2(Uout));     
        Ui2 = Hi.*Ui1;                
        Ui3 = ifft2(ifftshift(Ui2));     
        Uio = Ui3.*exp(1i.*screens(:,:,i));  
        I = Uio.*conj(Uio);              
        I_array(:,:,i) = I;            
        Uout = Uio;  
end

i = 2

i = 3

i = 4

i = 5

i = 6

i = 7

i = 8

i = 9

i = 10

i = 11

i = 12

i = 13

i = 14

i = 15

i = 16

i = 17

i = 18

i = 19

i = 20

i = 21

i = 22

i = 23

i = 24

i = 25

i = 26

i = 27

i = 28

i = 29

i = 30

i = 31

i = 32

i = 33

i = 34

i = 35

i = 36

i = 37

i = 38

i = 39

i = 40

i = 41

i = 42

i = 43

i = 44

i = 45

i = 46

i = 47

i = 48

i = 49

i = 50

i = 51

Answer 1

[Steven Lord]

Each time through the loop, you're completely overwriting the variable Uout. Compare:

for k = 1:5
    z = k.^2 % Overwrite with the new value
end

z = 1

z = 4

z = 9

z = 16

z = 25

with:

for k = 1:5
    x(k) = k.^2 % Add new value to the end
end

x = 1

x = 1×2

1 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x = 1×3

1 4 9

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x = 1×4

1 4 9 16

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x = 1×5

1 4 9 16 25

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Of course, since you know how large you want Uout to be, you should preallocate it before the loop starts so you fill it in rather than making it grow each iteration.

y = zeros(1, 5);
for k = 1:5
    y(k) = k.^2
end

y = 1×5

1 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 1×5

1 4 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 1×5

1 4 9 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 1×5

1 4 9 16 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y = 1×5

1 4 9 16 25

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The variable z is always a scalar since you're overwriting it each iteration.

The variable x changes its size as you add a new value to it at each iteration.

The variable y always has the same size even as each iteration changes one of its elements.

FYI the normalize function may be of interest to you in your construction of R.