I want to write a code to implement the the beam propagation algorithm for gaussian beams

I want the code to look at what happens when a Gaussian beam waist of 5mm passes through an element with a phase function exp(ikr, propagates for 100mm before passing through a 100mm focal length lens where an image is formed 100mm after the lens.

I want to output the intensity and phase of the beam various points during the propagation such as

(i) at the point just after the element (ii) the point just before the lens (iii) the point in the focal plane of the lens.

I want everything to be plotted in the x-z plane.

I know roughly what would happen however I am stuck at where to even start with writing the code for it.

Any help would be appreciated.

Here is what I've to atm

s = 4;
x=(-2*10^-s):(10^-(s+2)):(2*10^-s);   % x range

% d = 100*10^-3; % z distance after lens
f = 100*10^-3; % focus of the lens 

%ABCD matrix of thin lens with translation d 
A = 1;
B = f;
C = -1/f;
D = 0;

z = 100*10^-3; % propagation distance from source to lens 
lambda = 1*10^-6; % wavelength
k = 2*pi/lambda;
w0 = 5*10^-3; % beam waist

zR = (pi*(w0^2))/lambda; % Rayleigh range
w = w0*((1+(z/zR)^2)^.5); % spot size
R = z*(1+(zR/z)^2); % radius of curvature of the beam's wavefront

q1 = 1/((1/R)-((1i*lambda)/(pi*(w^2)))); % complex beam parameter before lens
q2 = (A*q1+B)/(C*q1+D); % complex beam parameter after lens

E = exp(-(((1i*k)/(2*q2)))*x.^2); % electric field
I = abs((E.^2))/2; % intensity 

figure(1)
plot(x,I)

[Image Analyst]

Do you have a mathematical equation for it?

[Jonathan]

No I don't, I thought about maybe using ABCD matrixes to map it but I dont know how to implement it in matlab