function morse_sofft
tic

% Position Grid Parameters for MO
xmin=-7;
xmax=7;
nx=256;
dx=(xmax-xmin)/(nx-1);
xgrid=linspace(xmin,xmax,nx);

% Momentum Grid Parameters for MO
p_xmin=-pi/dx;
p_xmax=-p_xmin;
p_xgrid=linspace(p_xmin,p_xmax-2*p_xmax/nx,nx);

% Potential Parameters for MO
mass=1.0;
De=150; %SC00:10.25
alpha=0.288; %SC00:0.2209
omega=alpha*sqrt(De*2/mass);

% Time step and initial Parameters for the Wave Function
delta_t=0.01;
steps=floor(10/delta_t)+1;
gamma_x=4.9883; %SC00:1
q_i=3.5; %SC00:5
p_i=0;
psi_x=cell(steps,1);
psi_p=cell(steps,1);
for l=1:steps
    psi_x{l}=zeros(nx,1);
    psi_p{l}=zeros(nx,1);
end
psi_x{1}=(gamma_x/pi)^0.25*exp(-0.5*gamma_x*(xgrid-q_i).^2+1i*p_i*(xgrid-q_i));

% plot(xgrid,abs(psi_x{1}));
psi_p{1}=fftshift(fft(psi_x{1}));
% define time grid for plotting purposes
time = zeros(steps,1);
% setup potential energy operator
potgrid=omega^2*xgrid.^2/2;%harmonic case
%potgrid=De*(1-exp(-alpha*xgrid)).^2;%Morse case
%potgrid=0.25*xgrid.^4;%quartic case
potential_operator=exp( -1i*potgrid*delta_t );
% setup kinetic energy (operator)
kingrid=0.5*(p_xgrid.^2/mass);
% kingrid( kingrid>50 ) = 50;
kingridshift=ifftshift( kingrid );
kinetic_e_operator=exp( -1i*kingridshift*delta_t );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% propagation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for l=2:steps
   
    time(l) = delta_t*(l-1);
    psi_x{l} = psi_x{l-1}.*potential_operator;
    psi_p{l} = fft(psi_x{l});
    psi_p{l} = psi_p{l}.*kinetic_e_operator;
    psi_x{l} = ifft(psi_p{l});
    
end
min(min((abs(psi_x{l})).^2));
max(max((abs(psi_x{l})).^2));

% for l=1:steps
%     plot(xgrid,(abs(psi_x{l})).^2); xlim([-5,10]); ylim([1.6013e-15,2.5]); title(['time: ', num2str(time(l))]); pause(0.01);
% end

auto(1)=1
x1(1)=q_i
for l=1:steps
    auto(l) = sum(psi_x{1}.*conj(psi_x{l}))*dx;
    x1(l)= sum(abs(psi_x{l}).^2.*xgrid)*dx;
    x2(l)= sum(abs(psi_x{l}).^2.*xgrid.^2)*dx;
    var(l)=x2(l)-x1(l)^2;
end
 
%output
 
plot(time,abs(auto));%autocorrelation
%plot(time,x1);%position
%plot(time,sqrt(var));%variance
toc
end