function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, V, t_idx, Observ)
    set(0,'defaulttextInterpreter','latex')
    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
    
    dt              =-1j*abs(this.TimeStepSize); % Imaginary Time
            
    KEop            = 0.5*(Transf.KX.^2+Transf.KY.^2);
    
    Observ.res      = 1; 
    AdaptIdx        = 0;

    % Change in Energy
    E               = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V);
    E               = E/Params.N;
    Observ.EVec     = [Observ.EVec E];

    % Chemical Potential
    muchem          = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
    Observ.mucVec   = [Observ.mucVec muchem];
                
    % Normalized residuals
    res             = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem);
    Observ.residual = [Observ.residual res];
    
    g_eff           = Params.gs      * (1/(sqrt(2*pi)*VParams.ell));
    gamma_eff       = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
    Ez              = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);

    while t_idx < Params.sim_time_cut_off

        % kin
        psi    = fftn(psi);
        psi    = psi.*exp(-0.5*1i*dt*KEop);
        psi    = ifftn(psi);
        
        % DDI
        frho   = fftn(abs(psi).^2);
        Phi    = real(ifftn(frho.*VDk));
        
        % Real-space
        psi    = psi.*exp(-1i*dt*(V + Ez + g_eff*abs(psi).^2 + gamma_eff*abs(psi).^3 + Params.gdd*Phi/(sqrt(2*pi)*VParams.ell) - muchem));
    
        % kin
        psi    = fftn(psi);
        psi    = psi.*exp(-0.5*1i*dt*KEop);
        psi    = ifftn(psi);
       
        % Renorm
        Norm   = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy;
        psi    = sqrt(Params.N)*psi/sqrt(Norm);
    
        muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
    
        if mod(t_idx,1000) == 0        
    
            % Change in Energy
            E               = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V);
            E               = E/Params.N;
            Observ.EVec     = [Observ.EVec E];
    
            % Chemical potential
            Observ.mucVec   = [Observ.mucVec muchem];
    
            % Normalized residuals
            res             = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem);
            Observ.residual = [Observ.residual res];
            Observ.res_idx  = Observ.res_idx + 1;
    
            %Adaptive time step -- Careful, this can quickly get out of control
            relres          = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);

            if relres <1e-4
                if AdaptIdx > 3 && abs(dt) > Params.mindt
                    dt = dt / 2;
                    AdaptIdx = 0;
                elseif AdaptIdx > 3 && abs(dt) < Params.mindt
                    break
                elseif Observ.residual(Observ.res_idx) < Params.rtol
                    break
                else
                    AdaptIdx = AdaptIdx + 1;
                end
            else
                AdaptIdx = 0;
            end
        end
        if any(isnan(psi(:)))
            disp('NaNs encountered!')
            break
        end
        t_idx=t_idx+1;
    end
    
    % Change in Energy
    E               = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V);
    E               = E/Norm;
    Observ.EVec     = [Observ.EVec E];
    
    % Normalized residuals
    res             = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem);
    Observ.residual = [Observ.residual res];
    Observ.res_idx  = Observ.res_idx + 1;
    
    % Chemical potential
    Observ.mucVec   = [Observ.mucVec muchem];
end