  In this work we review the theory of the spherical collapse model and
critically analyse the aspects of the numerical implementation of its
fundamental equations. By extending a recent work by Herrera et al. 2017, we
show how different aspects, such as the initial integration time, the
definition of constant infinity and the criterion for the extrapolation method
(how close the inverse of the overdensity has to be to zero at the collapse
time) can lead to an erroneous estimation (a few per mill error which
translates to a few percent in the mass function) of the key quantity in the
spherical collapse model: the linear critical overdensity $\delta_{\rm c}$,
which plays a crucial role for the mass function of halos. We provide a better
recipe to adopt in designing a code suitable to a generic smooth dark energy
model and we compare our numerical results with analytic predictions for the
EdS and the $\Lambda$CDM models. We further discuss the evolution of
$\delta_{\rm c}$ for selected classes of dark energy models as a general test
of the robustness of our implementation. We finally outline which modifications
need to be taken into account to extend the code to more general classes of
models, such as clustering dark energy models and non-minimally coupled models.
