function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
    
    % == User-defined potential == %
    V = this.Potential;
    assert(~anynan(V), 'Potential not defined! Specify as <SimulatorObject>.Potential = <PotentialsObject>.trap() + <AdditionalTerms>.');
    
    % == Calculating the DDIs == %
    if isfile(strcat(this.SaveDirectory, '/VDk_M.mat'))
        VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')));
        VDk = VDk.VDk;
    else
        VDk = this.Calculator.calculateVDkWithCutoff(Params,Transf,TransfRad);
        save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
    end
    fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType)
    
    % == Setting up the initial wavefunction == %
    psi = this.setupWavefunction(Params,Transf);
    
    if this.RunOnGPU
        psi = gpuArray(psi);
    end
end