Hello, dear forum members,
I have a question that I've been struggling with for some time, and despite extensive searching, I haven't been able to find an answer or identify the error.
I’m working with the discrete Fourier transform to solve the one- and two-dimensional Schrödinger equation using the split-operator method. This involves splitting the Hamiltonian into kinetic and potential parts. As long as the time steps are sufficiently small, the split-operator method is generally a suitable approach for determining the discrete wave function. In the one-dimensional case, this method works very well compared to others (e.g., Crank-Nicolson, see attached file). However, I’m encountering issues in the two-dimensional case, where I'm not obtaining the expected approximate wave function. I believe the error lies in the transformation to momentum space or I may have overlooked something basic.
In my example, a free wave packet with momentum k in the y-direction splits into two parts after a few time steps, which doesn't seem correct.
Could someone please help me with this issue? The example is taken from: J. M. Feagin, Quantum Methods with Mathematica, New York, TELOS/Springer-Verlag, 1994.
Thank you in advance!

TDFourier2D[timesteps_, potheight_, ellip_] := 
 Module[{NP, L, dx, dt, kx, ky, x0, y0, xmax, xmin, ymax, ymin, kxmax,
    kymax, \[Rho], TimeEnd, psi, K3, nt},
  NP = 100; (* netpoints *);
  L = 1;
  dx = L/(NP - 1); dt = 1/20 dx^2;
   TimeEnd = timesteps;
   \[Rho] = 0.16;  (* initial width *)
  x0 = 0.0 ; (* initial position *)
  y0 = 0.0;  (* initial position *)
  kx = 0 ;(* wave number x *)
  ky = 50; (* wave number y *)
  xmax = 1; (* dimension of the box *)
  xmin = -1; (* dimension of the box *)
  ymax = 1; (* dimension of the box *)
  ymin = -1; (* dimension of the box *)
  kxmax = N[Pi/dx]; (* momentum lattice x-direction*);
  kymax = N[Pi/dx];
  psi3[0] = 
   Array[1/(\[Rho] Sqrt[Pi]) Exp[-I kx #1 - (#1 - x0)^2/(2 \[Rho]^2)]*
       Exp[-I ky #2 - (#2 - y0)^2/(2 \[Rho]^2)] &, {NP, 
      NP}, {{xmin, xmax}, {ymin, ymax}}] // Quiet;
  K3 = Array[
    1/2 #1^2 + 1/2 #2^2 &, {NP, 
     NP}, {{-kxmax, kxmax}, {-kymax, kymax}}]; 
  V3 = Array[
    potheight/2 #1^2 + ellip*potheight/2 #2^2 &, {NP, 
     NP}, {{xmin, xmax}, {ymin, ymax}}]; 
  UV[dt] = Exp[-I V3 dt];
  UK[dt] = Exp[-I K3 dt];
  U[dt_]@psi_ := UV[dt]*InverseFourier[UK[dt]*Fourier[psi]];
  psi3[nt_, psi0_] := Nest[U[dt][#] &, psi0, nt];
  Do[psi3[nt] = psi3[1, psi3[nt - 1]], {nt, 1, TimeEnd, 1}];
  ]
ellipticalform = 
  0.2; (*elliptical form from from > 0 and <= 1; 1=cylindrical*)
TimeEnd = 160; (* timesteps 10 to 150 *)
potheight = 0; (*potential height, free particle *)
(TDFourier2D[TimeEnd, potheight, ellipticalform]) // AbsoluteTiming
graphTable1 = 
 GraphicsGrid[
  Partition[
   Table[Labeled[
     ListDensityPlot[Abs[psi3[i]], Axes -> True, PlotRange -> All, 
      ColorFunction -> "SunsetColors", 
      PerformanceGoal -> "Quality"], {"x", "y"}, {Left, Top}, 
     RotateLabel -> True], {i, 1, TimeEnd, TimeEnd/4}], 2], 
  Frame -> All, ImageSize -> 400]