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*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu));
    gamma_eff       = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2);
    EVar            = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu);

    pb = Helper.ProgressBar();
    pb.run('Running simulation in imaginary time: ');

    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 + EVar + g_eff*abs(psi).^2 + gamma_eff*abs(psi).^3 + Params.gdd*Phi/(sqrt(2*pi)*VParams.ell_eff) - 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;
    
            save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
            
            %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;
                    fprintf('Time step changed to '); disp(dt);
                    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];

    pb.run(' - Job Completed!');
    disp('Saving data...');
    save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
    disp('Save complete!');
end