  When $0<\alpha<1$, the approximation for the Caputo derivative
$$y^{(\alpha)}(x) = \frac{1}{\Gamma(2-\alpha)h^\alpha}\sum_{k=0}^n
\sigma_k^{(\alpha)} y(x-kh)+O\bigl(h^{2-\alpha}\bigr),$$ where
$\sigma_0^{(\alpha)} = 1, \sigma_n^{(\alpha)} = (n-1)^{1-a}-n^{1-a}$ and
$$\sigma_k^{(\alpha)} = (k-1)^{1-\alpha}-2k^{1-a}+(k+1)^{1-\alpha},\quad
(k=1...,n-1),$$ has accuracy $O\bigl(h^{2-\alpha}\bigr)$. We use the expansion
of $\sum_{k=0}^n k^\alpha$ to determine an approximation for the fractional
integral of order $2-\alpha$ and the second order approximation for the Caputo
derivative $$y^{(\alpha)}(x) = \frac{1}{\Gamma(2-\alpha)h^\alpha}\sum_{k=0}^n
\delta_k^{(\alpha)} y(x-kh)+O\bigl(h^{2}\bigr),$$ where $\delta_k^{(\alpha)} =
\sigma_k^{(\alpha)}$ for $2\leq k\leq n$, $$\delta_0^{(\alpha)} =
\sigma_0^{(\alpha)}-\zeta(\alpha-1), \delta_1^{(\alpha)} =
\sigma_1^{(\alpha)}+2\zeta(\alpha-1),\delta_2^{(\alpha)} =
\sigma_2^{(\alpha)}-\zeta(\alpha-1),$$ and $\zeta(s)$ is the Riemann zeta
function. The numerical solutions of the fractional relaxation and subdiffusion
equations are computed.
