  The time fractional ODEs are equivalent to convolutional Volterra integral
equations with completely monotone kernels. We therefore introduce the concept
of complete monotonicity-preserving ($\mathcal{CM}$-preserving) numerical
methods for fractional ODEs, in which the discrete convolutional kernels
inherit the $\mathcal{CM}$ property as the continuous equations. We prove that
$\mathcal{CM}$-preserving schemes are at least $A(\pi/2)$ stable and can
preserve the monotonicity of solutions to scalar nonlinear autonomous
fractional ODEs, both of which are novel. Significantly, by improving a result
of Li and Liu (Quart. Appl. Math., 76(1):189-198, 2018), we show that the
$\mathcal{L}$1 scheme is $\mathcal{CM}$-preserving, so that the $\mathcal{L}$1
scheme is at least $A(\pi/2)$ stable, which is an improvement on stability
analysis for $\mathcal{L}$1 scheme given in Jin, Lazarov and Zhou (IMA J.
Numer. Analy. 36:197-221, 2016). The good signs of the coefficients for such
class of schemes ensure the discrete fractional comparison principles, and
allow us to establish the convergence in a unified framework when applied to
time fractional sub-diffusion equations and fractional ODEs.
  The main tools in the analysis are a characterization of convolution inverses
for completely monotone sequences and a characterization of completely monotone
sequences using Pick functions due to Liu and Pego (Trans. Amer. Math. Soc.
368(12):8499-8518, 2016). The results for fractional ODEs are extended to
$\mathcal{CM}$-preserving numerical methods for Volterra integral equations
with general completely monotone kernels. Numerical examples are presented to
illustrate the main theoretical results.
