function fft_split_prop_henon_heiles_super_faster(hbar,lambda,x1,p1,x2,p2,zeta1,zeta2,M,dt,tmax,skipstep,wdisp)
%close all
%clear all
%M =256; %select even numbers
%if((zeta1==0.0)&&(zeta2==0.0))
L=16; %see the Henon-Heiles potential plot %range goes from -L/2 to +L/2
%else
%L=16;
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha1=(x1+1i*p1)/sqrt(2*hbar);
alpha2=(x2+1i*p2)/sqrt(2*hbar);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%making the state vectors%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = linspace(-L/2+L/M,L/2,M);
y=x;
if (zeta1==0.0)
psi_x=(1/(pi*hbar).^0.25)*exp(0.5*(alpha1*alpha1-abs(alpha1).^2))...
*exp(-(x/sqrt(2*hbar)-alpha1).^2);
else
psi_x=moller_squeeze_coherent_wavefunction_hbar(x1,p1,zeta1,x,hbar);
end
%I1= trapz(x,psi_x.*conj(psi_x))%%%
if (zeta2==0.0)
psi_y=(1/(pi*hbar).^0.25)*exp(0.5*(alpha2*alpha2-abs(alpha2).^2))...
*exp(-(y/sqrt(2*hbar)-alpha2).^2);
else
psi_y=moller_squeeze_coherent_wavefunction_hbar(x2,p2,zeta2,y,hbar);
end
psi_xy=psi_y.'*psi_x;
I12(1)= trapz(y,trapz(x,psi_xy.*conj(psi_xy)))%%%
%%%%%%%%%%%%%%%%Now the state vectors are created%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%Defining the Potential Grid and V(x,y)%%%%%%%%%%%%%%%%%%%%
[xx yy] = meshgrid(x,y);
Vxy=0.5*(xx.^2+yy.^2)+lambda*(xx.^2.*yy-yy.^3/3);
%surf(xx,yy,Vxy)
clear xx yy
%%%%%%%%%% Wavefunction is Transformed to momentum space%%%%%%%%%%%%%%%%%%%
psi_k = fft2(psi_xy);
%%%%%%%%%%%%%%%%%%%Fourrier coeficents%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
four_coef = [0:M/2 (-M/2+1):-1]; %fourrier coeficents
[coefx,coefy]=meshgrid(four_coef.^2,four_coef.^2);
halfcoefMat=exp((-1i*dt*hbar/4)*((2*pi/L).^2)*(coefx+coefy)); %half kinetic operation
revhalfcoefMat=exp((+1i*dt*hbar/4)*((2*pi/L).^2)*(coefx+coefy));%reverse half kinetic operation
fullcoefMat=exp((-1i*dt*hbar/2)*((2*pi/L).^2)*(coefx+coefy)); %full kinetic operation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
JJ=0;
%ttm=0.0;
%%%%%%%%%%%%%%taking exp(-idt/4*T)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
psi_k = psi_k.*halfcoefMat; %dt/2 propagator in momentum space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for indxt=0:tmax/dt
%%%%%%%%%%%%%%%%%%%%coming back to position space%%%%%%%%%%%%%%%%%%%%%%%%%%
psi_xy = ifft2(psi_k);
%%%%%%%%%%%%%%%%%%%%%multiply with the expoenential of the potential%%%%%%%
psi_xy=exp(-1i*(dt/hbar)*Vxy).*psi_xy;
%%%%%%%%%Now go to the momentum space to apply the exp(-idt/2T)%%%%%%%%%%%%
psi_k= fft2(psi_xy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
psi_k = psi_k.*fullcoefMat; %propagator in momentum space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(mod(indxt,skipstep)==0.0)
JJ=JJ+1;
%time(JJ)=indxt*dt;
%tmax/dt-indxt
psi_t=ifft2(psi_k.*revhalfcoefMat);
if (wdisp~=0)
pcolor(x,y,psi_t.*conj(psi_t)) %this is important to debug the program
%axis([-15 15 -15 15])
axis([-L/2 L/2 -L/2 L/2])
shading flat
drawnow()
end
end
end
%psi_xy = ifft2(psi_k);
%movie(MM)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end