This paper considers systems of M parabolic equations, of the form u dot i = S ijL(u j) , where L is an N-dimensional elliptic operator and the symmetric part of the coupling matrix ( S ij ) is positive semi-definite. A numerical scheme is explored in which each solution component is sequentially resolved according to an alternating direction implicit technique, requiring the solution of MN tridiagonal matrix equations per time-step. It is shown that unconditional stability is possible provided the asymmetries in the coupling matrix are not too large. In the case of purely symmetric coupling the stability conditions are shown to reduce to those for the single equation u ̇ = L(u) . Two applications of the method are discussed.