differential equations - Animated Wave propagation using split operator method

Here is the state of the art: By solving the time-dependent Schrödinger equation, one obtains that the propagated wave function after time step $\Delta t$ can be calculated by applying the time-evolution operator on the wave packet at any instant $t$ \begin{equation} \Psi({\vec r},t+\Delta t)=e^{-i H\Delta t}\Psi({\vec r},t), \end{equation} where we assumed that the Hamiltonian $H$ is time-independent.

Following the split-operator method, the time evolution operator for the Hamiltonian can be calculated as a series of matrix multiplications \begin{equation} \Psi({\vec r},t+\Delta t)=M1.M2.M1\Psi({\vec r},t)+O(\Delta t^3), \end{equation}

I would like to do iteratively:I start with psi(r,0) to calculate psi(r,Δt); then I use the psi(r,Δt), multiply it by the matrices and obtain psi(r,Δt+Δt); then I take psi(r,Δt+Δt), multiply it by the matrices and obtain psi(r,Δt+Δt+Δt), and so on.

 ass = {Δt>0, v > 0, b > 0, \hbar > 0, kx >= 0,ky >= 0};

 η[ℏ_, b_, x_, y_] =b/(2 Sqrt[π]) (x + I*y)*E^(-b*(x^2 + y^2)/(4*ℏ) )*{{1}, {1}};

M1[Δt_,v_,b_,y_,ℏ_] = {{1, 0}, {0, 1}}*Cos[(Δt*v*b*y)/(2 ℏ)]- 
      Sin[(Δt*v*b*y)/(2 ℏ)]/((Δt*v*b*y)/(2 ℏ))*{{0, -b y}, {-b y, 0}};



 M2[Δt_,v_,b_,y_,ℏ_,kx_,ky_] = {{1, 0}, {0, 1}}*Cos[(Δt*v Sqrt[kx^2 + ky^2])/(2 ℏ)]-Sin[Δt*v Sqrt[kx^2 + ky^2]]/(Δt*v Sqrt[kx^2 + ky^2])*{{0, kx - I*ky}, {kx + I*ky, 0}};

The Time Evolution of a Wave Packet can be calculated as a series of matrix multiplication

Ψ[Δt_,v_,b_,x_,y_,ℏ_, kx_,ky_] := M1[Δt, v, b,y, ℏ].M2[Δt, v, b, y, ℏ, kx,ky].M1[Δt, v, b, y, ℏ].η[ℏ, b,x, y];

This is the same expression as in Eq. (16) here