  The accuracy of the numerical solution of a fractional differential equation
depends on the differentiability class of the solution. The derivatives of the
solutions of fractional differential equations often have a singularity at the
initial point, which may result in a lower accuracy of the numerical solutions.
We propose a method for improving the accuracy of the numerical solutions of
the fractional relaxation and subdiffusion equations based on the fractional
Taylor polynomials of the solution at the initial point.
