function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)
% 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
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);
    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