Chemical Potential Calculator. See Landau and Lifshits, Stat.Phys. p. 163

n = 
        (10^9) ;   
g = 1  ;      
m = 23 (10^(-24)) ;  
k =  1.38 (10^(-16 ));  
hbar = 1.05 (10^(-27)) ;

Plot[μ /. 
        FindRoot[
          n == ((g (((m k T))^(3/2)))/((2^(1/
                          2)) (π^2) (hbar^3))) (∫_0^∞)((Sqrt[ z])/((E^(z - μ/(k T))) - 
                        1)) \[DifferentialD]z, {μ, -(10^(-15))}], {T, 1, 4}]

(* 
      Here we find the temperature where μ becomes  = 
        0 *) 
FindRoot[
      n == ((g (((m k T))^(3/2)))/((2^(1/
                      2)) (π^2) (hbar^3))) (∫_0^∞)((Sqrt[ z])/((E^z) - 1)) \[DifferentialD]z, {T, 
        (10^(-9))}]

(* 
      Or we can just use the formula *) 
T0 = (3.31/(g^(2/3))) ((hbar^2)/m) (n^(2/3))/
          k        (* L & 
            L (62.2) *)