diff --git a/GRIN_lens.m b/GRIN_lens.m
new file mode 100644
index 0000000..1ee7910
--- /dev/null
+++ b/GRIN_lens.m
@@ -0,0 +1,51 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; % wavelength in um
+N = 2^10; % number of points
+Lz = 600; % z-limit of visualization in um
+Lx = 600; % x-limit of visualization in um
+r = 50; % radius of the beam in um
+
+l_start = 150; % z-coordinate of the input surface of the lens
+l_end = 250; % z-coordinate of the output surface of the lens
+p = 0.01; % parameter of the lens
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+E = exp(-(x / r).^2); % Gaussian beam
+% E = rect(x / (2 * r)); % uniform beam
+
+l_width = l_end - l_start; % width of the gradient lens
+dn = @(x, z) -0.5 * p^2 * x.^2 * rect((z - l_start - (l_width / 2)) / l_width); % function of the GRIN lens
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dn(x, dz * n) * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end
+
+figure
+imagesc(z, x, rot90(I))
+colorbar
+title(sprintf("Beam radius = %.1f {\\mu}m. Lens start = %.1f {\\mu}m. Lens end = %.1f {\\mu}m. Lens parameter = %.3f", r, l_start, l_end, p))
+xlabel("z, mm")
+ylabel("x, mm")
+
+line([l_start l_start], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+line([l_end l_end], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+
+
+
+
+% colormap(gray)
\ No newline at end of file
diff --git a/axicon.m b/axicon.m
new file mode 100644
index 0000000..75a29da
--- /dev/null
+++ b/axicon.m
@@ -0,0 +1,45 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; %% wavelength in um
+N = 2^10; % number of points
+Lz = 200; % z-limit of visualization in um
+Lx = 600; % x-limit of visualization in um
+r = 100; % radius of the beam in um
+n0 = 1.5; % refractive index of the axicon
+a = 140; % apex angle of the axicon in degrees
+%%%%%
+
+k0 = 2*pi/lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+%% Axicon parameters %%
+g = 90 - a / 2;
+b = asind(n0 * cosd(a / 2)) + a / 2 - 90;
+p = sind(b);
+%%%%%
+
+tau = exp(1i * k0 * p * abs(x)); % phase function of the axicon
+E = exp(-x.^2 / r^2) .* tau; % Gaussian beam
+% E = rect(x / (2 * r)) .* tau; % uniform beam
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E).^2;
+end
+
+imagesc(z, x, rot90(I));
+title(sprintf("Beam radius = %.1f {\\mu}m. Apex angle = %.1f{\\deg}. Refractive index = %.1f", r, a, n0))
+xlabel("z, mm")
+ylabel("x, mm")
+colorbar
+%colormap(gray)
diff --git a/axicon_telescope.m b/axicon_telescope.m
new file mode 100644
index 0000000..51e1335
--- /dev/null
+++ b/axicon_telescope.m
@@ -0,0 +1,76 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1.064e-3; % wavelength in mm
+f1 = 80; % focal length of the 1st lens in mm
+f2 = 80; % focal length of the 2nd lens in mm
+N = 2^10; % number of points
+Lz = 2 * (f1 + f2); % z-limit of visualization in mm
+Lx = 10; %  x-limit of visualization in mm
+
+r = 1.5; % radius of the beam in mm
+n0 = 1.5; % refractive index of the axicon
+a = 176; % axicon apex angle in degrees
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+%% Axicon parameters %%
+g = 90 - a/2;
+b = asind(n0 * cosd(a / 2)) + a / 2 - 90;
+p = sind(b);
+%%%%%
+
+tau = exp(1i * k0 * p * abs(x)); % phase function of the axicon
+tau1 = exp(1i * k0 * x.^2 * 0.5 / f1); % phase function of the 1st lens
+tau2 = exp(1i * k0 * x.^2 * 0.5 / f2); % phase function of the 2nd lens
+
+E = exp(-(x / r).^2) .* tau; % Gaussian beam
+%E = rect(x / (2 * r)) .* tau; % uniform beam
+
+n = 1; % counter
+t = 0; % current z-coordinate
+
+I = zeros(N, N); % resulting intensity
+while t <= f1
+    E = exp(-1i * k0 * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E).^2;
+    t = t + dz;
+    n = n + 1;
+end
+
+E = E .* tau1;
+
+while t <= 2 * f1 + f2
+    E = exp(-1i * k0 * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E).^2;
+    t = t + dz;
+    n = n + 1;
+end
+
+E = E.*tau2;
+while t <= 2 * (f1 + f2)
+    E = ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E).^2;
+    t = t + dz;
+    n = n + 1;
+end
+
+imagesc(z, x, rot90(I), [0 4]);
+colorbar
+s_1 = sprintf("Kepler telescope with {f_1 =} %.2f mm and {f_2} = %.2f mm", f1, f2);
+s_2 = sprintf("Beam radius = %.1f mm. Apex angle = %.1f{\\deg}. Refractive index = %.1f", r, a, n0);
+title({s_1; s_2})
+xlabel("z, mm")
+ylabel("x, mm")
+line([f1 f1], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+line([2*f1+f2 2*f1+f2], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
diff --git a/bragg.m b/bragg.m
new file mode 100644
index 0000000..e37408c
--- /dev/null
+++ b/bragg.m
@@ -0,0 +1,49 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; % wavelength in um
+N = 2^9; % number of points
+Lz = 3000; % z-limit of visualization in um
+Lx = 1500; % x-limit of visualization in um
+r = 240; % radius of the beam in um
+
+phi = 7 % incident angle of the beam in degrees
+T = lambda / (2 * sind(phi)); % period of the grating in um
+% T = 5 * lambda;
+
+g_start = 500; % z-coordinate of the input surface of the grating
+g_end = 1000; % z-coordinate of the output surface of the grating
+%%%%%
+k0 = 2*pi/lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+tau = exp(-1i * k0 * sind(phi) * x); % phase function of the grating
+E = exp(-(x / r).^2) .* tau;
+%E = rect(x / (2*r)) .* tau; %uniform beam
+
+g_width = g_end - g_start; % width of the gradient grating
+dn = @(x, z) 0.001 * rect((z - g_start - (g_width / 2)) / g_width) * cos(2 * pi / T * x);
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dn(x, dz*n) * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end   
+
+figure
+imagesc(z, x, rot90(I))
+colorbar
+title(sprintf("Beam radius = %.1f {\\mu}m. Incident angle = %.1f{\\deg} Grating start = %.1f {\\mu}m. Grating end = %.1f {\\mu}m. Grating period = %.3f {\\mu}m", r, phi, g_start, g_end, T))
+xlabel("z, mm")
+ylabel("x, mm")
+
+line([g_start g_start], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+line([g_end g_end], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
\ No newline at end of file
diff --git a/gradient_glass.m b/gradient_glass.m
new file mode 100644
index 0000000..303a1ae
--- /dev/null
+++ b/gradient_glass.m
@@ -0,0 +1,48 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; % wavelength in um
+N = 2^10; % number of points
+Lz = 2000; % z-limit of visualization in um
+Lx = 600; % x-limit of visualization in um
+r = 60; % radius of the beam in um
+
+g_start = 500; % z-coordinate of the input surface of the glass
+g_end = 1500; % z-coordinate of the output surface of the glass
+p = 0.07; % parameter of the glass
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+E = exp(-(x / r).^2); % Gaussian beam
+%E = rect(x / (2 * r)); % uniform beam
+
+g_width = g_end - g_start; % width of the gradient glass
+dn = @(x, z) p * rect((z - g_start - (g_width / 2)) / g_width) * (x / Lx); % function of the gradient glass
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dn(x, dz*n) * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end   
+
+figure
+imagesc(z, x, rot90(I))
+colorbar
+title(sprintf("Beam radius = %.1f {\\mu}m. Glass start = %.1f {\\mu}m. Glass end = %.1f {\\mu}m. Medium parameter = %.3f", r, g_start, g_end, p))
+xlabel("z, mm")
+ylabel("x, mm")
+
+line([g_start g_start], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+line([g_end g_end], [-Lx/2 Lx/2], 'Color','red','LineStyle','--');
+
+% colormap(gray)
\ No newline at end of file
diff --git a/kerr_effect.m b/kerr_effect.m
new file mode 100644
index 0000000..724eec6
--- /dev/null
+++ b/kerr_effect.m
@@ -0,0 +1,38 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; % wavelentgh in um
+N = 2^9; % number of points
+Lz = 3000; % z-limit of visualization in um
+Lx = 300; % x-limit of visualization in um
+A = 0.8; % amplitude of the beam
+r = 30; % radius of the beam in um;
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+E = A * exp(-(x / r).^2); % Gaussian beam
+% E = A * rect(x / (2 * r)); % uniform beam
+
+dn = @(E) (10^-3) * abs(E).^2; % change of refractive index vs field
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dn(E) * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end
+
+imagesc(z, x, rot90(I));
+title(sprintf("beam amplitude = %.1f. Beam radius = %.1f mm", A, r))
+xlabel("z, mm")
+ylabel("x, mm")
+colorbar
\ No newline at end of file
diff --git a/lens.m b/lens.m
new file mode 100644
index 0000000..b62d3af
--- /dev/null
+++ b/lens.m
@@ -0,0 +1,41 @@
+clc
+clear
+
+%% User choices %%
+x_max = 2; % % maximum x in mm
+r = 2; % radius of the beam in mm
+f = 100; % focal length of the lens
+lambda = 1e-3; % wavelemgth in mm
+N = 2^10; % number of points
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+Lz = 2 * f; % z-limit of visualization in mm
+Lx = 2 * x_max; % x-limit of visualization in mm
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+tau = exp(1i * k0 * x.^2 * 0.5 / f); % lens phase function
+
+E = exp(-(x / r).^2) .* tau; % Gaussian beam
+% E = rect(x / (2 * r)) .* tau; % Uniform beam
+
+I = zeros(N, N); % resulting intensity
+
+for n = 1:N
+    E = ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end
+
+imagesc(z, x, rot90(I));
+title(sprintf("Focal length = %.1f mm. Beam radius = %.1f mm", f, r))
+xlabel("z, mm")
+ylabel("x, mm")
+colorbar
+% colormap(gray)
diff --git a/rect.m b/rect.m
new file mode 100644
index 0000000..56f8ff8
--- /dev/null
+++ b/rect.m
@@ -0,0 +1,3 @@
+function r = rect(x);
+  r = 0.5 * (sign(x + 0.5) - sign(x - 0.5));
+end
\ No newline at end of file
diff --git a/telescope.m b/telescope.m
new file mode 100644
index 0000000..2e5a819
--- /dev/null
+++ b/telescope.m
@@ -0,0 +1,69 @@
+clc
+clear
+
+%% User choices %%
+x_max = 2; % maximum x in mm
+f1 = 100; % focal length of the 1st lens of the telescope in mm
+f2 = 100; % focal length of the 2nd lens of the telescope in mm
+r = 0.5; % radius of the beam in mm
+lambda = 1e-3; % wavelength in mm
+N = 2^10; % number of points
+%%%%%
+
+k0 = 2 * pi / lambda; % wavenumber
+Lz = 2 * (f1 + f2); % z-limit of visualization in mm
+Lx = 2 * x_max; % x-limit of visualization in mm
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz; % z coordinate array
+x = -Lx/2:dx:Lx/2; % x coordinate array
+kx = -pi/dx:dkx:pi/dx; % kx coordinate array
+
+tau1 = exp(1i * k0 * x.^2 * 0.5 / f1); % 1st lens phase function
+tau2 = exp(-1i * k0 * x.^2 * 0.5 / f2); % 2nd lens phase function
+
+E = exp(-(x / r).^2); % Gaussian beam
+% E = rect(x / (2 * r)); % Uniform beam
+
+t = 0; % iterator for z
+n = 1; % counter
+I = zeros(N, N); % resulting intensity 
+
+while t <= f1
+    E = exp(-1i * k0 * dz) * ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E); 
+    t = t + dz;
+    n = n + 1;
+end
+
+E = E .* tau1;
+
+while t <= 2 * f1 + f2
+    E = exp(-1i * k0 * dz) * ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E); 
+    t = t + dz;
+    n = n + 1;
+end
+
+E = E.*tau1;
+
+while t <= 2 * (f1 + f2)
+    E = exp(-1i * k0 * dz) * ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E); 
+    t = t + dz;
+    n = n + 1;
+end
+
+figure
+imagesc(z, x, rot90(I))
+line([f1 f1], [-x_max x_max], 'Color','red','LineStyle','--');
+line([2*f1+f2 2*f1+f2], [-x_max x_max], 'Color','red','LineStyle','--');
+
+colorbar
+title(sprintf("Kepler telescope with {f_1 =} %.2f mm and {f_2} = %.2f mm", f1, f2))
+xlabel("z, mm")
+ylabel("x, mm")
+% colormap(gray)
diff --git a/two_beams.m b/two_beams.m
new file mode 100644
index 0000000..47876c2
--- /dev/null
+++ b/two_beams.m
@@ -0,0 +1,52 @@
+clc
+clear
+
+%% User choices %%
+lambda = 1; % wavelength in um
+k0 = 2 * pi / lambda; % wavenumber
+N = 2^10; % number of points
+Lz = 1000; % z-limit of visualization in um
+Lx = 1200; % x-limit of visualization in um
+
+x0_1 = 200; % x-position of the 1st beam in um
+x0_2 = -200; % x-position of the 2nd beam in um
+
+r_1 = 60; % radius of the 1st beam in um
+r_2 = 60; % radius of the 2nd beam in um
+
+t_1 = 10; % tilt of the 1st beam in degrees
+t_2 = 10; % tilt of the 2nd beam in degrees
+%%%%%
+
+dz = Lz / (N - 1);
+dx = Lx / (N - 1);
+dkx = 2 * pi / Lx;
+
+z = 0:dz:Lz;
+x = -Lx/2:dx:Lx/2;
+kx = -pi/dx:dkx:pi/dx;
+
+%% Gaussian beam %%
+E1 = exp(-((x - x0_1) / r_1).^2 + 1i * k0 * sind(t_1) * x); 
+E2 = exp(-((x - x0_2) / r_2).^2 - 1i * k0 * sind(t_2) * x);
+%% Uniform beam %%
+% E1 = rect((x - x0_1) / (2 * r_1)) .* exp(1i * k0 * sind(t_1) * x); 
+% E2 = rect((x - x0_2) / (2 * r_2)) .* exp(-1i * k0 * sind(t_2) * x);
+%%%%%
+E = E1 + E2;
+
+I = zeros(N, N); % resulting intensity
+for n = 1:N
+    E = exp(-1i * k0 * dz) .* ifft(fftshift(exp(-1i * sqrt(k0^2 - kx.^2) * dz) .* fftshift(fft(E))));
+    I(n, :) = abs(E);    
+end   
+
+imagesc(z, x, rot90(I));
+s_1 = "Two crossing beams";
+s_2 = sprintf("First beam: {x_0} = %.1f {\\mu}m; radius = %.1f {\\mu}m; tilt = %.1f{\\deg}", x0_1, r_1, t_1);
+s_3 = sprintf("Second beam: {x_0} = %.1f {\\mu}m; radius = %.1f {\\mu}m; tilt = %.1f{\\deg}", x0_2, r_2, t_2);
+title({s_1; s_2; s_3})
+xlabel("z, {\\mu}m")
+ylabel("x, {\\mu}m")
+colorbar
+% colormap(gray)
\ No newline at end of file