The study of the spectrum of coupled random matrices has received rather little attention. To the best of our knowledge, coupled random matrices have been studied, to some extent, by Mehta. In this work, we explain how the integrable technology can be brought to bear to gain insight into the nature of the distribution of the spectrum of coupled Hermitean random matrices and the equations the associated probabilities satisfy. In particular, the two-Toda lattice, its algebra of symmetries and its vertex operators will play a prominent role in this interaction. Namely, the method is to introduce time parameters, in an artificial way, and to dress up a certain matrix integral with a vertex integral operator, for which we find Virasoro-like differential equations. These methods lead to very simple nonlinear third-order partial differential equations for the joint statistics of the spectra of two coupled Gaussian random matrices.Comment: 56 pages, published version, abstract added in migratio