Abstract. We discuss the theoretical convergence and numerical evaluation of simultaneous interpolation quadrature formulas which are exact for rational functions. Basically, the problem consists in integrating a single function with respect to different measures by using a common set of quadrature nodes. Given a multi-index n, the nodes of the integration rule are the zeros of the multiorthogonal Hermite-Padé polynomial with respect to (S, α, n), where S is a collection of measures, and α is a polynomial which modifies S. The theory is based on the connection between Gausslike simultaneous quadrature formulas of rational type and multipoint Hermite-Padé approximation. As for the numerical treatment we present a procedure based on the technique of modifying the integrand by means of a change of variable when the integrand has real poles close to the integration interval. The output of some tests show the power of this approach when compared to other ones.Key words. simultaneous integration rule, Gauss-like rational quadrature formula, meromorphic integrand AMS subject classifications. Primary 41A55. Secondary 41A28, 65D321. Introduction. Simultaneous integration is a problem which arises in computer graphics to determinate the color of light which emanates from a given point on a surface toward the viewer (see Borges [4] and the bibliography therein). In an abstract setting, this problem was earlier studied by Nikishin [19] in connection with Hermite-Padé approximation. Simultaneous integration means that we are integrating a single function f : I → R, with respect to m distinct measures ds 1 , ..., ds m on I, respectively.A reference index Φ = (Ov + 1)/N u (performance ratio) is considered in [4] to state the efficiency of the procedure. Here Ov is the overall degree of exactness and N u is the number of integrand evaluations. Borges [4] remarked that when m ≥ 3 the use of the m corresponding Gaussian rules of polynomial type yields Φ < 1, which indicates a low performance. Instead he suggests quadrature rules whose nodes are the zeros of multi-orthogonal polynomials for which Φ > 1 holds.Geometric rate of convergence has been proved by the authors in [8] for polynomials methods when the integrand is analytic in a neighborhood of I. Nevertheless, instability shows up when the corresponding numerical method is applied to meromorphic integrands with poles close to I. Several methods to integrate functions with singularities on or near the integration interval have been developed in the past few years.An interesting approach, whose starting point seems to be [13], is that based on the use of rational Gaussian integration rules connected with multipoint Padé