We consider the shift-invariant space, S(Phi), generated by a set Phi = {phi(1),...,phi(r)} of compactly supported distributions on R when the vector of distributions phi := (phi(1),..., phi(r))(T) satisfies a system of refinement equations expressed in matrix form as phi = Sigma(alpha is an element of Z) a(alpha)phi(2 . - alpha) where a is a finitely supported sequence of r x r matrices of complex numbers. Such multiple refinable functions occur naturally in the study of multiple wavelets. The purpose of the present paper is to characterize the accuracy of Phi, the order of the polynomial space contained in S(Phi), strictly in terms of the refinement mask a. The accuracy determines the L-p-approximation order of S(Phi) when the functions in Phi belong to L-p(R) (see Jia [10]). The characterization is achieved in terms of the eigenvalues and eigenvectors of the subdivision operator associated with the mask a. In particular, they extend and improve the results of Heil, Strang and Strela [7], and of Plonka [16]. In addition, a counterexample is given to the statement of Strang and Strela [20] that the eigenvalues of the subdivision operator determine the accuracy. The results do not require the linear independence of the shifts of phi.