Résumé: In this paper, we develop a new numerical method to estimate a multivariate probit model. To this end, we derive a new decomposition of normal multivariate integrals that has two appealing properties. First, the decomposition may be written as the sum of normal multivariate integrals, in which the highest dimension of the integrands is reduced relative to the initial problem. Second, the domains of integration are bounded and delimited by the correlation coefficients. Application of a Gauss-Legendre quadrature rule to the exact likelihood function of lower dimension allows for a major reduction of computing time while simultaneously obtaining consistent and efficient estimates for both the slope and the scale parameters. A Monte Carlo study shows that the finite sample and asymptotic properties of our method compare extremely favorably to the maximum simulated likelihood estimator in terms of both bias and root mean squared error. [Authors]