function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
    
    format long;
    
    % Convergence Criteria:
    epsilon         = 1E-8;
    alpha           = 1E-4;
    beta            = 0.9;

    Observ.residual = 1; 
    Observ.res      = 1;
    psi_old         = psi;         % Previous psi value (for heavy-ball method)
    
    if this.PlotLive
        Plotter.plotLive(psi,Params,Transf,Observ)
        drawnow
    end

    % Minimization Loop
    for idx      = 1:this.MaxIterationsForGD
        
        % Compute gradient
        J        = compute_gradient(psi, Params, Transf, VDk, V);

        % Calculate chemical potential and norm
        muchem   = sum(real(conj(psi(:)) .* J(:))) / sum(abs(psi(:)).^2);
        
        % Calculate residual and check convergence
        residual = sum(abs(J(:) - muchem * psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
        if residual < epsilon
            fprintf('Convergence reached at iteration %d\n', idx);
            break;
        end
        
        % Update psi using heavy-ball method
        psi_new  = (1 + beta) * psi - alpha * J - beta * psi_old;
        psi_old  = psi;
        
        % Normalize psi
        Norm     = sum(abs(psi_new(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
        psi      = sqrt(Params.N) * psi_new / sqrt(Norm);

        if mod(idx,100) == 0        
    
            % Collect change in energy
            E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
            E               = E/Norm;
            Observ.EVec     = [Observ.EVec E];
    
            % Collect chemical potentials
            Observ.mucVec   = [Observ.mucVec muchem];
    
            % Collect residuals
            Observ.residual = [Observ.residual residual];
            Observ.res_idx  = Observ.res_idx + 1;
            
            if this.PlotLive
                Plotter.plotLive(psi,Params,Transf,Observ)
                drawnow
            end

            save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
        end
    end
    
    % Check if max iterations were hit without convergence
    last_iteration_number = idx;
    if last_iteration_number == this.MaxIterationsForGD
        fprintf('Max iterations reached without convergence. Final chemical potential: %.6f, Residual: %.6f\n', muchem, residual);
    else
        fprintf('Converged in %d iterations. Final chemical potential: %.6f\n', last_iteration_number, muchem);
    end
    
    % Change in Energy
    E              = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
    E              = E/Norm;
    Observ.EVec    = [Observ.EVec E];

    disp('Saving data...');
    save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
    disp('Save complete!');
    
end

%% Helper functions
% Numerical Gradient Calculation using the finite differences method
function J      = compute_gradient(psi, Params, Transf, VDk, V)

    % Operators
    
    % Kinetic energy
    KEop    = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
    HKin    = @(w) real(ifft(KEop.*fft(w)));

    % Trap Potential
    HV      = @(w) V.*w;

    % Contact interactions
    Hint    = @(w) (Params.gs*abs(psi).^2).*w;
    
    % DDIs
    frho    = fftn(abs(psi).^2);
    Phi     = real(ifftn(frho.*VDk));
    Hddi    = @(w) (Params.gdd*Phi).*w;

    % Quantum fluctuations
    Hqf     = @(w) (Params.gammaQF*abs(psi).^3).*w;
    
    H       = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);

    J = H(psi);

end