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.
