We develop an efficient numerical scheme for the 3D mean-field spherical dynamo equation. The scheme is based on a semi-implicit discretization in time and a spectral method in space based on the divergence-free spherical harmonic functions. A special semi-implicit approach is proposed such that at each time step one only needs to solve a linear system with constant coefficients. Then, using expansion in divergence-free spherical harmonic functions in the transverse directions allows us to reduce the linear system at each time step to a sequence of one-dimensional equations in the radial direction, which can then be efficiently solved by using a spectral-element method. We show that the solution of fully discretized scheme remains bounded independent of the number of unknowns, and present numerical results to validate our scheme.