We develop a multidimensional Stein methodology for non-degenerate self-decomposable random vectors in $\mathbb{R}^d$ having finite first moment. Building on previous univariate findings, we solve an integro-partial differential Stein equation by a mixture of semigroup and Fourier analytic methods. Then, under a second moment assumption, we introduce a notion of Stein kernel and an associated Stein discrepancy specifically designed for infinitely divisible distributions. Combining these new tools, we obtain quantitative bounds on smooth-Wasserstein distances between a probability measure in $\mathbb{R}^d$ and a non-degenerate self-decomposable target law with finite second moment. Finally, under an appropriate spectral gap assumption, we investigate, via variational methods, the existence of Stein kernels. In particular, this leads to quantitative versions of classical results on characterizations of probability distributions by variational functionals.