Summary: Multiple Hermite polynomials are an extension of the classical Hermite polynomials for which orthogonality conditions are imposed with respect to \(r>1\) normal (Gaussian) weights \(w_j(x)=e^{-x^2+c_jx}\) with different means \(c_j/2\), \(1 \leq j \leq r\). These polynomials have a number of properties such as a Rodrigues formula, recurrence relations (connecting polynomials with nearest neighbor multi-indices), a differential equation, etc. The asymptotic distribution of the (scaled) zeros is investigated, and an interesting new feature is observed: depending on the distance between the \(c_j\), \(1 \leq j \leq r\), the zeros may accumulate on \(s\) disjoint intervals, where \(1 \leq s \leq r\). We will use the zeros of these multiple Hermite polynomials to approximate integrals of the form \(\int_{-\infty}^{\infty} f(x) \exp(-x^2 + c_jx)\, dx\) simultaneously for \(1 \leq j \leq r\) for the case \(r=3\) and the situation when the zeros accumulate on three disjoint intervals. We also give some properties of the corresponding quadrature weights.