Bi-orthogonal polynomials and related objects (quasi-polynomials and their Fourier-Laplace transforms) play an essential role in determining the correlation functions in the two-matrix model.  In the case of a certain class of measures defined by two polynomial potentials $V_1$ and $V_2$ one can deduce sequences of systems of first order finite-dimensional ODEs for these bi-orthogonal polynomials with polynomial coefficients.  In contrast to the case of orthogonal polynomials, where one finds 2-dimensional system of ODEs, here the analogous objects consist of two pairs of ODEs, each pair of size equal to the degree of the potentials $V_1$, $V_2$

The associated spectral curves of these four systems are all equal (a spectral "quadrality"), a novel feature which is absent in the 1-matrix case where there is a single relevant ODE.

Each of the four systems admits a Lax-pair representation of the discrete string equations and compatible isomonodromy deformation equations.

In fact we can associate to any of the four system a natural Riemann-Hilbert problem which is technically complicated by the fact that the leading matrix coefficient at the irregular singularity $x = \infty$ has a highly degenerate spectrum in which all but one eigenvalue vanish.

The precise formulation of the RH problem will be the key to opening the way to the analysis of the large $n$ asymptotics by means of the inverse monodromy method already successfully applied by A. Its and P. Bleher in the case of one-matrix models.

In collaboration with B. Eynard and J. Harnad.