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

Question

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.

Answer 1

If the different numerical solutions do not converge upon refinement, then there is an error of some sort. If I were in your position I would check the numerical solutions against some known analytical solution. There should be some solutions to the linear SE if I remember correctly. I would start with those and see if your methods already disagree for the simple case and then work your way up in complexity. Also, is there some global physical properties like total energy/impulse etc. that you can check?