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 (Numer. Math. 67, 271–288 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 (J. Comput. Appl. Math. 178, 131–145 2005) but it seems to have gone unnoticed. We describe the result in detail for \(r=2\) and give some examples.
