The generalized Fermi-Dirac functions and their derivatives are important in evaluating the thermodynamic quantities of partially degenerate electrons in hot dense stellar plasmas. New recursion relations of the generalized Fermi-Dirac functions have been found. An effective numerical method to evaluate the derivatives of the generalized Fermi-Dirac functions up to third order with respect to both degeneracy and temperature is then proposed, following Aparicio. A Fortran program based on this method, together with a sample test case, is provided. Accuracy and domain of reliability of some other, popularly used analytic approximations of the generalized Fermi-Dirac functions for extreme conditions are investigated and compared with our results.