using GLMakie
nbose(x) = x > 0 ? 1 / (exp(x) - 1) : NaN
nfermi(x) = 1 / (exp(x) + 1)
nmb(x) = exp(-x)
ϵ = range(-5, 5, 1000)

fig = Figure(size = (800, 400),fontsize=28,linewidth=4,fonts= (; regular = "CMU Serif",bold="CMU Serif"),)
ax = Axis(fig[1, 1], xlabel =L"(\epsilon-\mu)/k_BT",ylabel = L"\langle n\rangle",
xticks=LinearTicks(5),
yticks=[0,1,2],)

lines!(ax, ϵ, nbose.(ϵ), color = :red,label="Bose-Einstein")
lines!(ax, ϵ, nmb.(ϵ), color = :mediumseagreen,label="Maxwell-Boltzmann")
lines!(ax, ϵ, nfermi.(ϵ), color = :dodgerblue,label="Fermi-Dirac")
ylims!(ax,0, 2.5)
xlims!(ax,extrema(ϵ)...)
axislegend(ax,framevisible=false)
hidedecorations!(ax,label=false,ticks=false,ticklabels=false)
fig