diff --git a/Quasi2DBogoliubovSpectrum.m b/Quasi2DBogoliubovSpectrum.m
new file mode 100644
index 0000000..2d96620
--- /dev/null
+++ b/Quasi2DBogoliubovSpectrum.m
@@ -0,0 +1,75 @@
+%% Physical constants
+PlanckConstant               = 6.62607015E-34;
+PlanckConstantReduced        = 6.62607015E-34/(2*pi);
+FineStructureConstant        = 7.2973525698E-3;
+ElectronMass                 = 9.10938291E-31;
+GravitationalConstant        = 6.67384E-11;
+ProtonMass                   = 1.672621777E-27;
+AtomicMassUnit               = 1.660539066E-27; 
+BohrRadius                   = 5.2917721067E-11; 
+BohrMagneton                 = 9.274009994E-24; 
+BoltzmannConstant            = 1.38064852E-23;
+StandardGravityAcceleration  = 9.80665;
+SpeedOfLight                 = 299792458;
+StefanBoltzmannConstant      = 5.670373E-8;
+ElectronCharge               = 1.602176634E-19;
+VacuumPermeability           = 1.25663706212E-6; 
+DielectricConstant           = 8.8541878128E-12;
+ElectronGyromagneticFactor   = -2.00231930436153;
+AvogadroConstant             = 6.02214076E23;
+ZeroKelvin                   = 273.15;
+GravitationalAcceleration    = 9.80553;
+VacuumPermittivity           = 1 / (SpeedOfLight^2 * VacuumPermeability);
+HartreeEnergy                = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
+AtomicUnitOfPolarizability   = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3
+
+% Dy specific constants
+Dy164Mass                    = 163.929174751*1.660539066E-27;
+Dy164IsotopicAbundance       = 0.2826;
+DyMagneticMoment             = 9.93*9.274009994E-24;
+
+%% Dispersion relation of the quasiparticle excitations
+AtomNumber        = 1E5;
+wz                = 2*pi*72.4;
+lz                = sqrt(PlanckConstantReduced/(Dy164Mass*wz));                                      % Defining a harmonic oscillator length
+as                = 102.4*BohrRadius;                                                                % Scattering length
+Trapsize          = 7.6;
+alpha             = 0;
+phi               = 0;
+MeanWidth         = 2.8215042184E3*lz;
+k                 = linspace(0, 1e7, 1000);
+
+
+AtomNumberDensity = AtomNumber / (Trapsize * lz)^2;
+add               = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
+eps_dd            = add/as;
+gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
+
+[fk,Fka,Ukk] = computePotentialInMomentumSpace(k, lz, alpha, phi, gs, eps_dd);
+
+% == Quantum Fluctuations term == %
+gQF               = ((256 * PlanckConstantReduced^2) / (15*Dy164Mass*MeanWidth^3)) * as^(5/2) * (1 + ((3/2) * eps_dd^2));
+
+
+DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
+EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 
+
+figure(1)
+set(gcf,'Position',[100 100 950 750])
+% xvals = (k .* lz/sqrt(2));
+xvals = (k .* add);
+yvals = EpsilonK ./ (PlanckConstantReduced * wz);
+plot(xvals, yvals,LineWidth=2.0)
+% xlim([3.45, 3.65])
+% ylim([0, 0.001])
+title(horzcat(['$a_s = ',num2str(1/eps_dd),'a_{dd}, '], ['na_{dd}^2 = ',num2str(AtomNumberDensity * add^2),'$']),'fontsize',16,'interpreter','latex')
+xlabel('$ka_{dd}$','fontsize',16,'interpreter','latex')
+ylabel('$\epsilon(k)/\hbar \omega_z$','fontsize',16,'interpreter','latex')
+grid on
+
+%%
+function [fk,Fka,Ukk] = computePotentialInMomentumSpace(k, lz, alpha, phi, gs, eps_dd)
+    fk                = (3 * sqrt(pi)) * (k .* lz/sqrt(2)) .* exp((k .* lz/sqrt(2)).^2) .* erfc((k .* lz/sqrt(2))) ;
+    Fka               = (fk .* sin(deg2rad(phi))^2 - 1) + (cos(deg2rad(alpha))^2 .* (3 - (fk .* (sin(deg2rad(phi))^2 + 1))));
+    Ukk               = (gs/ (sqrt(2 * pi) * lz)) .* (1 + (eps_dd .* Fka));
+end
\ No newline at end of file