We present an algorithm for the numerical solution of (possibly nonlinear) fractional differential equations of the form y(α) (t) = f (t,y(t), y(β1) (t), y(β2) (t), ..., y(βn) (t)) with α > βn > βn-1 > ⋯ > β1 > 0, combined with suitable initial conditions. The fractional derivatives are described in the Caputo sense. The algorithm is based on Adomian's decomposition approach and the solutions are calculated in the form of a convergent series with easily computable components. Several numerical examples are presented. © 2006 Elsevier Inc. All rights reserved.