In this chapter, an adaptation of Stein’s method for bounding the error in multivariate normal approximation is presented. For simplicity, the distance measure used is based on expectations of functions with three bounded derivatives; more natural measures of distance would require much more complicated treatments. The Stein equation used is now a second-order partial differential equation. Solutions to the equation are exhibited, and some of their properties are established; they can then be used to derive a general bound on the approximation error in multivariate normal approximation. For exploiting the general bound, a local approach is introduced, which uses a multivariate version of the double decomposition used for (univariate) normal approximation. This is applied to the number of monochrome edges in a graph whose vertices are randomly coloured. A size bias coupling approach is also developed and applied to the joint distribution of counts of vertices of different degrees in the Bernoulli random graph.