Hi Babr, I think this is basically what is intended, taking the plot out of the for loop.
close all; clear; clc
N = 2048;
L = 1e4;
dx = 4e-3;
[x0,y0] = meshgrid((-N/2:N/2-1)*dx);
[xL,yL] = meshgrid((-N/2:N/2-1)*dx);
[fx,fy] = meshgrid((-N/2:N/2-1)/(N*dx));
Ni = 36;
xi= 1e-2 .*[-62.5,-62.5,-62.5,-62.5,-62.5,-62.5,-37.5,-37.5,-37.5,-37.5,-37.5,-37.5,...
-12.5,-12.5,-12.5,-12.5,-12.5,-12.5,12.5,12.5,12.5,12.5,12.5,12.5,37.5,37.5,37.5,37.5,37.5,37.5,...
62.5,62.5,62.5,62.5,62.5,62.5,];
yi= 1e-2 .*[62.5,37.5,12.5,-12.5,-37.5,-62.5,62.5,37.5,12.5,-12.5,-37.5,-62.5,62.5,37.5,12.5,-12.5,-37.5,-62.5,...
62.5,37.5,12.5,-12.5,-37.5,-62.5,62.5,37.5,12.5,-12.5,-37.5,-62.5,62.5,37.5,12.5,-12.5,-37.5,-62.5,];
w0 = 10e-2;
P = 1e3;
lambda = 1080e-9; k = 2*pi/lambda;
L0= 10; l0= 5e-3;
kL= 2.*pi./L0; kl= 5.92./l0;
Cn2= randi(1,60)*1e-16;
Kx = 2.*pi.*fx; Ky = 2.*pi.*fy;
U0= 0*x0;
for idx=1:Ni
Utemp = sqrt(2.*P./pi)./w0 .* exp(-((x0-xi(idx)).^2+(y0-yi(idx)).^2)./w0.^2);
Utemp(((x0-xi(idx)).^2+(y0-yi(idx)).^2)>w0.^2) = 0;
U0 = U0 + Utemp;
end
figure(1)
s = pcolor(x0, y0, abs(U0));
s.EdgeColor = 'none';
title("0^{th} Intenesity screen")
xlabel 'X [m] '
ylabel 'Y [m] '
cbar = colorbar;
cbar.Label.String = 'Intensity(X, Y, 0) [W/cm^2]';
ax = gca;
axis square
ax.FontName = 'Times new roman';
ax.FontSize = 15;
