  Burchnall and Chaundy showed that if two ODOs $P$, $Q$ with analytic
coefficients commute there exists a polynomial $f(\lambda ,\mu)$ with complex
coefficients such that $f(P,Q)=0$, called the BC-polynomial. This polynomial
can be computed using the differential resultant for ODOs. In this work we
extend this result to matrix ordinary differential operators, MODOs. Matrices
have entries in a differential field $K$, whose field of constants $C$ is
algebraically closed and of zero characteristic. We restrict to the case of
order one operators $P$, with invertible leading coefficient. A new
differential elimination tool is defined, the matrix differential resultant. It
is used to compute the BC-polynomial $f$ of a pair of commuting MODOs and
proved to have constant coefficients. This resultant provides the necessary and
sufficient condition for the spectral problem $PY=\lambda Y \ , \ QY=\mu Y$ to
have a solution. Techniques from differential algebra and Picard-Vessiot theory
allow us to describe explicitly isomorphisms between commutative rings of MODOs
$C[P,Q]$ and a finite product of rings of irreducible algebraic curves.
