diff --git a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
new file mode 100644
index 0000000..db7ed45
--- /dev/null
+++ b/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
@@ -0,0 +1,48 @@
+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;
+
+gs_pf_2D     = 1/(sqrt(2*pi)*VParams.sigma);
+gdd_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));
+
+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);
+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;
+
+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);
+
+end
+