Burchnall and Chaundy showed that if two ODOs P , Q with analytic coefficients commute there exists a polynomial f (λ, µ) 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 P Y = λY , QY = µ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.