function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)

    gs           = Params.gs; 
    gdd          = Params.gdd;
    gammaQF      = Params.gammaQF; 
    
    KEop         = 0.5*(Transf.KX.^2+Transf.KY.^2);
    
    g_pf_2D      = 1/(sqrt(2*pi)*VParams.ell);
    gQF_pf_2D    = sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2));
    Ez           = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
    muchem_tilde = muchem - Ez;
    
    % eigs only works with column vectors
    psi          = psi.';
    KEop         = KEop.';
    VDk          = VDk.';
    
    % Interaction Potential
    frho         = fftn(abs(psi).^2);
    Phi          = real(ifftn(frho.*VDk));
    % Operators
    H            = @(w) real(ifft(KEop.*fft(w)));
    C            = @(w) (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF*abs(psi).^3.*w);
    muHC         = @(w) (-muchem_tilde * w) + H(w) + C(w);
    X            = @(w) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;

    % 2-D matrix will be unravelled to a single column vector and the corresponding BdG matrix of (N^2)^2 elements solved for.
    Size         = length(psi(:)); 
    Neigs        = length(psi(:));
    opts.tol     = 1e-16;
    opts.disp    = 1;
    opts.issym   = 0;
    opts.isreal  = 1;
    opts.maxit   = 1e4;
    
    BdGVec       = @(g) BdGSolver2D.BdGMatrix(g, psi, Params, VDk, VParams, Transf, muchem); % This function takes a column vector as input and returns a
                                                                                             % matrix-vector product which is also a column vector
    [g,D]        = eigs(BdGVec,Size,Neigs,'sr',opts);
    evals        = diag(D); 
    clear D;
    
    % Eigenvalues
    evals = sqrt(evals);

    % Obtain f from g
    for ii = 1:Neigs
        gres    = reshape(g(:,ii), size(psi));
        f(:,ii) = reshape((1/evals(ii)) * (muHC(gres) + (2.*X(gres))), [], 1);
    end
    
    % Obtain u and v from f and g
    u = (f + g)/2;
    v = (f - g)/2;
    
    % Renormalize to \int |u|^2 - |v|^2 = 1
    for ii=1:Neigs
        normalization = sum(abs(u(:,ii)).^2 - abs(v(:,ii)).^2);
        u(:,ii)       = u(:,ii)/sqrt(normalization);
        v(:,ii)       = v(:,ii)/sqrt(normalization);
    end
    
    modes.u = u'; modes.v = v';
    modes.g = g'; modes.f = f';
    
end