  In this article, we study the numerical approximation of stochastic
differential equations driven by a multidimensional fractional Brownian motion
(fBm) with Hurst parameter greater than 1/3. We introduce an implementable
scheme for these equations, which is based on a second order Taylor expansion,
where the usual Levy area terms are replaced by products of increments of the
driving fBm. The convergence of our scheme is shown by means of a combination
of rough paths techniques and error bounds for the discretisation of the Levy
area terms.
