A sparse-grid method for solving multi-dimensional backward stochastic differential equations (BSDEs) based on a multi-step time discretization scheme [31] is presented. In the multi-dimensional spatial domain, i.e. the Brownian space, the conditional mathematical expectations derived from the original equation are approximated using sparse-grid Gauss-Hermite quadrature rule and (adaptive) hierarchical sparse-grid interpolation. Error estimates are proved for the proposed fully-discrete scheme for multi-dimensional BSDEs with certain types of simplified generator functions. Finally, several numerical examples are provided to illustrate the accuracy and efficiency of our scheme.Mathematics subject classification: 60H10, 60H35, 65C10, 65C20, 65C50.Because analytical solutions of BSDEs are often very difficult to obtain, approximate numerical solutions of BSDEs become highly desired in practical applications. There are mainly two types of numerical methods for BSDEs. One is based on the relation between the forwardbackward stochastic differential equations (FBSDEs) and corresponding parabolic partial differential equation (PDEs) [13,14,22]; the other is directly based on BSDEs or FBSDEs [2,3,6,7,10,12,23,27,28,30,31]. Zhao et al. proposed a θ-scheme for BSDEs in [28]; in [29], it was extended to a generalized θ-scheme. In [31], a stable multi-step scheme was proposed which is a highly accurate numerical method for BSDEs. Note that for the second type of numerical methods, approximating spatial derivatives at different time-space points for the case of PDEs is converted to approximating conditional mathematical expectations with Gaussian kernels centered at different time-space points.It should be noted that the BSDEs used in practice usually involve a multi-dimensional Brownian motion, such as the option pricing problem with multiple underlying assets. Existing numerical methods for BSDEs can be theoretically extended to the multi-dimensional cases; however, the computational cost may be unaffordable due to the so-called curse of dimensionality. The most popular approach to solving multi-dimensional BSDEs is the Monte Carlo method [3,28] that is very easy to implement. However, the convergence rate is typically very slow, although having a mild dependence on the dimensionality. Thus, an accurate and efficient numerical method for solving multi-dimensional BSDEs is highly desired in the BSDE community.In this paper, we extend the multi-step method in [31] using the sparse-grid method for solving multi-dimensional BSDEs. As discussed in [31], the target BSDE (1.1) is discretized by the multi-step scheme in the time direction. In the spatial domain, a quadrature rule is needed to approximate all the conditional mathematical expectations (multi-dimensional integrals) and an interpolation scheme is also needed to evaluate the integrands of the expectations at nongrid quadrature points. The sparse-grid method is highly suitable for the multi-step scheme because it has been demonstrated to be effective and efficient in dea...