This paper is concerned with the numerical solution of the initial value problem for the fractional differential equation of order \( \beta, \) \(( 0 0 \), with \( u - u_0 = 0, t \leq 0\) where \( \Phi \) is a sufficiently regular function and \( D^{\beta}\) is the fractional differential operator in the sense of \textit{Caputo} [Geophys. J. R. Astr. Soc. 13, 529--539 (1967)]. To obtain a scheme that approximates the solution of this problem, the authors introduce the fractional integrator operator \( I^{\beta}\) of order \( \beta \) such that \( I^{\beta} \cdot D^{\beta} \) is the identity operator and they construct a mixed collocation method on \( u(t) - u_0 = I^{\beta} ( \Phi ( u (\cdot), t ) ) \) with a basis of ``hat'' functions similar to the used in mixed finite elements. Some numerical experiments, with linear and non linear \( \Phi \)-functions and fractional order \( \beta = 1/2, \) are presented. In particular, interesting numerical results are obtained for Joulin model of spherical flames that possess a singularity whose detection depends on the grid size.