A numerical algorithm for calculating the mean energy of radiotherapy electron beams has been developed. This algorithm is fast and accurate which makes it suitable for routine clinical use. First, a Gaussian distribution of the electron energy spectrum is derived from the linear Boltzmann equation. Based on this Gaussian spectrum, and after introducing a correction to the CSDA mean energy, a recursive-iterative algorithm for the mean energy calculation is developed. The multiple-scattering correction is taken into account using Yang's pathlength distribution theory. Numerical results of the present algorithm are compared with the results obtained through Monte Carlo simulation as well as Harder's formula. Good agreement with Monte Carlo simulation is achieved. Also the new algorithm is much more accurate than the commonly-used empirical formula of Harder.