A Golub-Welsch version for simultaneous Gaussian quadrature

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. 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.