In this paper, we investigate the asymptotic distribution of the normalized error for the Mittag--Leffler Euler (MLE) method applied to a class of multidimensional fractional stochastic differential equations. These equations are reformulated as stochastic Volterra equations (SVEs) featuring a non-diagonal, matrix-valued kernel $K(u)=u^{\alpha-1}E_{\alpha,\alpha}(Au^{\alpha})$ with singular exponent $\alpha \in (\frac{1}{2}, 1)$. To enhance computational efficiency, the singular kernel is discretized using the left-rectangle rule, posing technical challenges for the theoretical analysis. To address this, we introduce an auxiliary $K$-undiscretized scheme to bridge the gap between the exact solution and the MLE method, integrating Jacod's stable convergence theory for conditional Gaussian martingales with methodologies developed for SVEs. To the best of our knowledge, this is the first work to establish the asymptotic error distribution for numerical methods incorporating non-diagonal matrix-valued kernels.