  The solutions of fractional differential equations (FDEs) have a natural
singularity at the initial point. The accuracy of their numerical solutions is
lower than the accuracy of the numerical solutions of FDEs whose solutions are
differentiable functions. In the present paper we propose a method for
improving the accuracy of the numerical solutions of ordinary linear FDEs with
constant coefficients which uses the fractional Taylor polynomials of the
solutions. The numerical solutions of the two-term and three-term FDEs are
studied in the paper.
