diff --git a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
index 8458d8d..cc26b1e 100644
--- a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
+++ b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
@@ -5,7 +5,7 @@ gs           = Params.gs;
 gdd          = Params.gdd;
 gammaQF      = Params.gammaQF; 
 
-KEop         = 0.5*kx.^2 + 0.5*ky.^2;
+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;
 
@@ -23,15 +23,15 @@ 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;
+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
+% Operate in order on g
 BdG          = @(g) muHC(muHC(g) + (2.*X(g)));
 
 syssize      = size(psi);
-opts.v0      = psi(:);
+opts.v0      = psi(:,:);
 opts.tol     = 1e-16;
 opts.disp    = 1;
 opts.issym   = 0;
@@ -46,17 +46,10 @@ 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))));
+    f(:,:,ii) = (1/evals(ii)) * (muHC(g(:,:,ii)) + (2.*X(g(:,:,ii))));
 end
 
 % Obtain u and v from f and g
@@ -65,9 +58,9 @@ 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);
+    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';