In recent decades, numerous and varied numerical methods have been proposed and studied to approximate solutions for various classes of fractional differential equations, primarily those involving single-term or multiple-order equations. However, equations incorporating fractional iterated derivatives have not received widespread attention. In this work we describe a reliable strategy to approximate the solution of higher-order fractional differential equations where both the fractional derivative and the iterated derivatives are described in the Caputo sense. Specifically, we propose a fractional variational iteration method (FVIM) where the Lagrange multiplier associated with the correction term is explicitly determined by means of the Laplace transform.

For the second-order case, we give a sufficient condition -involving the coefficients of the equation and the fractional order of the Caputo derivative- which guarantees the convergence of the sequence generated by the FVIM. Furthermore, this convergence is independent of the initial function considered for the iteration.

Finally, some examples are presented in order to illustrate the applicability of the method and the reliability of the theoretical results obtained. In particular, for most of them we observe that the FVIM leads to the exact solution which shows the power of the method in practice.