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 ---
    VarArray                  = Params.VarArray;
    VParams.ell               = VarArray(1); 
    VParams.nu                = VarArray(2);
    [VParams]                 = this.Calculator.find_nk_fit(VParams);
    
    % --- Initialize --- 
    mkdir(sprintf(this.SaveDirectory))
    mkdir(sprintf('./Data/Run_%03i',Params.njob))
    fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8);
    

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

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

        VParams.ell           = VarArray(1);
        VParams.nu            = VarArray(2);
        [VParams]             = this.Calculator.find_nk_fit(VParams);
        
        [psi,V,VDk]           = this.initialize(Params,VParams,Transf);
        
        % --- Adding some noise ---
        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 + 4*noise;

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

        save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',runIdx,nn),'psi','Observ','Transf','Params','VDk','V','VarArray');

        if relelldiff < Params.ellcutoff && relnudiff < Params.nucutoff
            break
        end
    end
    save(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','Observ','Transf','Params','VDk','V','VarArray');
    delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',runIdx))
    
end