  We study the Crank-Nicolson scheme for stochastic differential equations
(SDEs) driven by multidimensional fractional Brownian motion $(B^{1}, \dots,
B^{m})$ with Hurst parameter $H \in (\frac 12,1)$. It is well-known that for
ordinary differential equations with proper conditions on the regularity of the
coefficients, the Crank-Nicolson scheme achieves a convergence rate of
$n^{-2}$, regardless of the dimension. In this paper we show that, due to the
interactions between the driving processes $ B^{1}, \dots, B^{m} $, the
corresponding Crank-Nicolson scheme for $m$-dimensional SDEs has a slower rate
than for the one-dimensional SDEs. Precisely, we shall prove that when $m=1$
and when the drift term is zero, the Crank-Nicolson scheme achieves the exact
convergence rate $n^{-2H}$, while in the case $m=1$ and the drift term is
non-zero, the exact rate turns out to be $n^{-\frac12 -H}$. In the general case
when $m>1$, the exact rate equals $n^{\frac12 -2H}$. In all these cases the
limiting distribution of the leading error is proved to satisfy some linear SDE
driven by Brownian motions independent of the given fractional Brownian
motions.
