  We present an exponentially convergent numerical method to approximate the
solution of the Cauchy problem for the inhomogeneous fractional differential
equation with an unbounded operator coefficient and Caputo fractional
derivative in time. The numerical method is based on the newly obtained
solution formula that consolidates the mild solution representations of
sub-parabolic, parabolic and sub-hyperbolic equations with sectorial operator
coefficient $A$ and non-zero initial data. The involved integral operators are
approximated using the sinc-quadrature formulas that are tailored to the
spectral parameters of $A$, fractional order $\alpha$ and the smoothness of the
first initial condition, as well as to the properties of the equation's
right-hand side $f(t)$. The resulting method possesses exponential convergence
for positive sectorial $A$, any finite $t$, including $t = 0$ and the whole
range $\alpha \in (0,2)$. It is suitable for a practically important case, when
no knowledge of $f(t)$ is available outside the considered interval $t \in [0,
T]$. The algorithm of the method is capable of multi-level parallelism. We
provide numerical examples that confirm the theoretical error estimates.

Numerical Analysis;Mathematical Software;Numerical Analysis;Analysis of PDEs;Classical Analysis and ODEs