Occupation number

(* occupation number at temperature T, chemical potential μ, energy ϵ; s = 1 (fermi), 0 (Boltzmann), or -1 (bose) *)
Clear[f];
f[T_,μ_,ϵ_,s_]:=(1/((E^((ϵ-μ)/T))+s));

(* fermi, Boltzmann, bose *)
Block[{T=1,μ=-1,s={1,0,-1}},Plot[Evaluate[f[T,μ,ϵ,s]],{ϵ,0,4},PlotRange->All]]

(* fermi *)
Block[{T=1,μ={-1,0,1,10},s=1},Plot[Evaluate[f[T,μ,ϵ,s]],{ϵ,0,20},PlotRange->All]]