The theory used to model the behavior of the beam is describded in P. A. Bélanger, "Beam propagation and the ABCD ray matrices," Opt. Lett. 16, 196-198 (1991).
The code computes each plane starting from the far left to the far right (the direction of propagation through the cylinder) and computes both direction independently (in calcul_x and calcul_y respectively). For each plane, it computes the field then the intensity. After getting the two intensity planes, a threshold (found experimentally for the sillica, 1/e^2 for the first version) is given to be used as the radii (semi-minor and semi-major axis) of the ellipse (in calcul_ellipse). Each ellipse is drawn to the 3D plot one after the other (again following the propagation direction).
For further details on the implementation, please refer to the code which is commented extensively.
You can find images (intensity and isometric view) from different angles (XZ plane, YZ planes) for different depths here.