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.