We consider the following Schrödinger system

$$\{\begin{array}{rl}& -\mathrm{\Delta }{u}_j=\sum _{i=1}^k{\beta }_{ij}|u_i{|}^3|u_j|u_j+{\lambda }_j|u_j{|}^{q-2}{u}_j,\text{in}\mathrm{\Omega },\\ & u_j=0\text{on}\mathrm{\partial }\mathrm{\Omega },j=1,\cdots ,k\end{array}$$

where $\mathrm{\Omega }\subset {\mathbb{R}}^3$ is a bounded domain with smooth boundary. Assume $5<q<6,{\lambda }_j>0,{\beta }_{jj}>0,j=1,\cdots ,k$, ${\beta }_{ij}={\beta }_{ji},i\ne j,i,j=1,\cdots ,k$. Note that the nonlinear coupling terms are of critical Sobolev growth in dimension $3$. We prove that under an additional condition on the coupling matrix the problem has infinitely many sign-changing solutions. The result is obtained by combining the method of invariant sets of descending flow with the approach of using approximation of systems of subcritical growth.