# Solving the pulse propagation using four different FDTD methods gives four different results - Which to trust?

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.