InvFD, an OCTAVE routine for the numerical inversion of the Fermi-Dirac integral

Given a value of the Fermi-Dirac integral $\displaystyle \int_0^{\infty}\frac{t^j\,dt}{e^{t-\eta}+1}$, this routine returns the value of the parameter η for a given j with a relative error less than 10^{-13}. The inversion is provided for η ∈ [-30,100] and for the following set of j values, { -1/2, 1/2, 3/2, 5/2, 1, 2 }. For η ∈ [-30,-2], an iterative method involving the McDougall and Stoner series is used. For the rest of the region of interest, a rational fit is employed.