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

wz_tilde     = Params.wz / Params.w0;
gs           = Params.gs; 
gdd          = Params.gdd;
gammaQF      = Params.gammaQF; 

KEop         = 0.5*kx.^2 + 0.5*ky.^2;
Ez           = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2);
muchem_tilde = muchem - Ez;

g_pf_2D      = 1/(sqrt(2*pi)*VParams.sigma);
gQF_pf_2D    = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2));

% 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) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF.* abs(psi).^3).*w;

% Operate on g
BdG          = @(g) muHC(muHC(g) + (2.*X(g)));

syssize      = size(psi);
opts.v0      = psi(:);
opts.tol     = 1e-16;
opts.disp    = 1;
opts.issym   = 0;
opts.isreal  = 1;
opts.maxit   = 1e4;
Neigs        = syssize;

[g,D]        = eigs(BdG,syssize,Neigs,'sr',opts);
evals        = diag(D); 
clear D;

% Eigenvalues
evals = sqrt(evals);

BogModeCutoff = 1;
if BogModeCutoff == 1
    g(:,real(evals)<1e-6) = [];
    evals(real(evals)<1e-6) = [];
    Neigs = length(evals);
end

% Obtain f from g
f = zeros(size(g));
for ii = 1:Neigs
    f(:,ii) = (1/evals(ii)) * (muHC(g(:,ii)) + (2.*X(g(:,ii))));
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)*Transf.dz;
    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