E = 1000;

b = pi/4:0.00001:pi/3;
nmax = floor(E./b - 0.5);
dos = 0*b;

for counter = 1:length(b),
    n = 0:nmax(counter);
    dos(counter) = b(counter)*sum(1./sqrt(E - b(counter)*(n + 0.5)));
    dos(counter) = dos(counter)/(2*sqrt(E));
end

figure(1)
plot(E./b,dos)
axis([1100 1120 0.96 1.06])
grid minor on
xlabel('E/(hbar omega_c)')
ylabel('DOS(B)/DOS(B=0)')

figure(2)
plot(E./b,dos)
axis([1000 1200 0.9 2])
xlabel('E/(hbar omega_c)')
ylabel('DOS(B)/DOS(B=0)')
