  The zeros of type II multiple orthogonal polynomials can be used for
quadrature formulas that approximate $r$ integrals of the same function $f$
with respect to $r$ measures $\mu_1,\ldots,\mu_r$ in the spirit of Gaussian
quadrature. This was first suggested by Borges in 1994, even though he does not
mention multiple orthogonality. We give a method to compute the quadrature
nodes and the quadrature weights which extends the Golub-Welsch approach using
the eigenvalues and left and right eigenvectors of a banded Hessenberg matrix.
This method was already described by Coussement and Van Assche in 2005 but it
seems to have gone unnoticed. We describe the result in detail for $r=2$ and
give some examples.
