We combine Stein's method with Malliavin calculus in order to obtain explicit bounds in the multidimensional normal approximation (in the Wasserstein distance) of functionals of Gaussian fields. Our results generalize and refine the main findings by Peccati and Tudor (2005), Nualart and Ortiz-Latorre (2007), Peccati (2007) and Nourdin and Peccati (2007b, 2008); in particular, they apply to approximations by means of Gaussian vectors with an arbitrary, positive definite covariance matrix. Among several examples, we provide an application to a functional version of the Breuer-Major CLT for fields subordinated to a fractional Brownian motion.