  The aim of this paper is to describe a Matlab package for computing the
simultaneous Gaussian quadrature rules associated with a variety of multiple
orthogonal polynomials. Multiple orthogonal polynomials can be considered as a
generalization of classical orthogonal polynomials, satisfying orthogonality
constraints with respect to $r$ different measures, with $r \ge 1$. Moreover,
they satisfy $(r+2)$--term recurrence relations. In this manuscript, without
loss of generality, $r$ is considered equal to $2$. The so-called simultaneous
Gaussian quadrature rules associated with multiple orthogonal polynomials can
be computed by solving a banded lower Hessenberg eigenvalue problem.
Unfortunately, computing the eigendecomposition of such a matrix turns out to
be strongly ill-conditioned and the \texttt{Matlab} function \texttt{balance.m}
does not improve the condition of the eigenvalue problem. Therefore, most
procedures for computing simultaneous Gaussian quadrature rules are implemented
with variable precision arithmetic. Here, we propose a \texttt{Matlab} package
that allows to reliably compute the simultaneous Gaussian quadrature rules in
floating point arithmetic. It makes use of a variant of a new balancing
procedure, recently developed by the authors of the present manuscript, that
drastically reduces the condition of the Hessenberg eigenvalue problem.
