I'd like to simulate the propagation of a pulse, and have different options for solving that.
On the one hand, I can use the non-linear schrödinger equation
$$\partial_zE=\frac{i}{2k_0}\nabla^2_\perp E$$
with $k_0$ the propagation constant. I can solve it either in the cartesian coordinate system (using a fourier transform) or in the cylindrical coordinate system (using the hankel transform).
On the other hand, I can use the non-paraxial pulse propagation equation, as defined in https://doi.org/10.1103/PhysRevLett.117.043902
$$\partial_z\hat{E}=ik_z\hat{E}$$
with $k_z$ defined either as a paraxial approximation in
$$k_z = k_0 - \frac{k_x^2+k_y^2}{2k_0}$$
or the non-paraxial version
$$k_z = \sqrt{k_0^2-k_x^2-k_y^2}$$
with $k_x$ and $k_y$ the spatial frequencies.
Finally, I also can calculate the pulse and it's behavior for a gaussian beam.
I did that, and plotted everything in one figure:
The non-paraxial propagation with the paraxial constant overlaps directly with the paraxial approximation in the cartesian system.
Still, I get four different results for five different methods. Which of them is correct, and which of them is not? A change of the step size or resolution in x/y-direction does not change the result in a significant way.
My final goal is to use the correct solution for automated tests later, but without knowing which solution I will get this is difficult.