We evaluate the fractional integral by making use of the integral representation of the convolution kernel. We construct an efficient quadrature of this integral representation and use it as a part of a fast time stepping method. The (possible) singularity of f is taken into account. This new method is particularly well-suited for long time simulations.