diff --git a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
index db7ed45..11e1757 100644
--- a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
+++ b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
@@ -9,8 +9,8 @@ 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;
 
-gs_pf_2D     = 1/(sqrt(2*pi)*VParams.sigma);
-gdd_pf_2D    = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2));
+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.';
@@ -21,11 +21,13 @@ VDk          = VDk.';
 frho         = fftn(abs(psi).^2);
 Phi          = real(ifftn(frho.*VDk));
 
+% Operators
 H            = @(w)     real(ifft(KEop.*fft(w)));
-C            = @(w)     (((gs_pf_2D*gs*abs(psi).^2) + (gdd_pf_2D*gdd*Phi)).*w) + (gammaQF.* abs(psi).^3 .*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)*(gammaQF.* abs(psi).^3).*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);
@@ -44,5 +46,32 @@ 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) = 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