  The well-known Caputo fractional derivative and the corresponding Caputo
fractional integral occur naturally in many equations that model physical
phenomena under inhomogeneous media. The relationship between the two
fractional terms can be readily obtained by applying the Laplace transform to a
given equation. We seek to numerically approximate Caputo fractional integrals
using a Taylor series expansion for convolution integrals. This naturally
extends into being able to approximate convolution integrals for a wider class
of convolution integral kernels $K(t-s)$. One of the main advantages under this
approach is the ability to numerically approximate weakly singular kernels,
which fail to converge under traditional quadrature methods. We provide
stability and convergence analysis for these composite quadratures, which offer
optimal convergence for approximating functions in $C^{\gamma}[0,T]$, where
$\alpha \leq \gamma \leq 5$ and $0<\alpha < 1$. For the order $\gamma =
1,2,3,4,5$ scheme, the resulting approximation is $O(\tau^{\gamma})$ accurate,
where $\tau$ is the size of the partition of the time domain. By instead
utilizing a fractional Taylor series expansion, we are able to obtain for
$\gamma \in (0,5)-\{1,2,3,4\}$ order scheme, which yields an approximation of
$O(\tau^{\gamma})$ with a constant dependent on the kernel function which
improves the order of convergence. This allows for a far wider class of
functions to be approximated, and by strengthening the regularity assumption,
we are able to obtain more accurate results. General convolution integrals
exhibit these results up to $\gamma = 2$ without the assumption of $K$ being
decreasing. Finally, some numerical examples are presented, which validate our
findings.
