We consider the shift-invariant space, $\bbbs(\Phi)$, generated by a set $\Phi=\{\phi_1,\ldots,\phi_r\}$ of compactly supported distributions on $\RR$ when the vector of distributions $\phi:=(\phi_1,\ldots,\phi_r)^T$ satisfies a system of refinement equations expressed in matrix form as $$ \phi=\sum_{\alpha\in\ZZ}a(\alpha)\phi(2\,\cdot - \,\alpha) $$ where $a$ is a finitely supported sequence of $r\times r$ matrices of complex numbers. Such {\it multiple refinable functions} occur naturally in the study of multiple wavelets. The purpose of the present paper is to characterize the {\it accuracy} of $\Phi$, the order of the polynomial space contained in $\bbbs(\Phi)$, strictly in terms of the refinement mask $a$. The accuracy determines the $L_p$-approximation order of $\bbbs(\Phi)$ when the functions in $\Phi$ belong to $L_p(\RR)$ (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$.