Burchnall and Chaundy showed that if two ordinary differential operators (ODOs) P, Q with analytic coefficients commute then 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. Our 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. We define a new differential elimination tool, the matrix differential resultant. We use it to compute the BC-polynomial f of a pair of commuting MODOs, and we also prove that it has constant coefficients. This resultant provides the necessary and sufficient condition for the spectral problem PY=λ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.