function [Params, Transf, psi, V, VDk] = run(this)
    format long;

    % --- Obtain simulation parameters --- 
    [Params]                  = this.setupParameters();
    this.SimulationParameters = Params;
    
    % --- Set up spatial grids and transforms --- 
    [Transf]                  = this.setupSpace(Params);

    % --- Set up for Variational method ---
    VParams.SelfConIter       = Params.SelfConIter;      % Max number of iterations to perform self-consistent calculation
    VParams.ell               = Params.ell;              % initial [ell], ell is the "width" - psi ~ e^(z^2/ell^2)
                                  
    % Window of optimization 
    VParams.ell_lower         = Params.ell_lower;
    VParams.ell_upper         = Params.ell_upper;
    
    % Relative cutoffs
    VParams.ellcutoff         = Params.ellcutoff;
    
    % --- Initialize --- 
    mkdir(sprintf(this.SaveDirectory))
    mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob))
    fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8);

    [psi,V,VDk]               = this.initialize(Params,VParams,Transf);
    ells(1)                   = VParams.ell; 
    E_Var                     = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
    E_vs_iter(1)              = E_Var(ells(1));

    for nn = 1:Params.SelfConIter
        Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; Observ.res_idx = 1;
        
        t_idx                 = 1; % Start at t = 0;

        [~,V,VDk]             = this.initialize(Params,VParams,Transf); % Do not overwrite psi, susbequent iterations should use psi generated at the end of the loop to converge faster
                                                                        % This is still needed however to recalculate VDk with new value
                                                                        % for the variational parameter
                                                                        
        % --- Adding some noise --- 
        % Noise added in every iteration to ensure it is not stuck in some local minimum
        Norm                  = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation
        psi                   = sqrt(Params.N)*psi/sqrt(Norm);
        r                     = normrnd(0,1,size(psi));
        theta                 = rand(size(psi));
        noise                 = r.*exp(2*pi*1i*theta);
        psi                   = psi + Params.nsf*noise; 

        % --- Run --- 
        [psi, Observ]         = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn);
        psi                   = gather(psi);
        
        % --- Constrained minimization --- 
        E_Var                 = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
        VParams.ell           = fmincon(E_Var,VParams.ell,[],[],[],[],Params.ell_lower,Params.ell_upper,[],fminconoptions);
        
        % --- Convergence check --- 
        ells                  = [ells VParams.ell]; 
        relelldiff            = abs(ells(nn+1)-ells(nn))/ells(nn); 
        E_vs_iter             = [E_vs_iter E_Var(VParams.ell)];

        save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs_%i.mat'),Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
        
        %Plotting
        if this.PlotLive
            figure(2);
            plot(ells,E_vs_iter,'LineStyle','none','Marker','o','MarkerSize',3,'Color','b','MarkerFaceColor','b')
            xlabel('$\ell$', 'FontSize', 14)
            ylabel('$E_{var}$', 'FontSize', 14)
            title('Variational Energy', 'FontSize', 14);
            drawnow
        end
        
        if relelldiff < Params.ellcutoff
            break
        end
    end
    fprintf('\n')
    disp('Saving data...');
    save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
    disp('Save complete!');
    delete(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs_*.mat'),Params.njob))

end