<abstract><p>We consider the following Schrödinger system</p>
<p><disp-formula> <label/> <tex-math id="FE1"> \begin{document}$ \begin{equation*} \left\{ \begin{aligned} &-\Delta u_j = \sum\limits_{i = 1}^k \beta_{ij}|u_i|^3|u_j|u_j+\lambda_j|u_j|^{q-2}u_j, \ \ \ \text{in}\, \, \Omega, \\ &u_j = 0\quad\text{on}\, \, \partial\Omega, \, \, j = 1, \cdots, k \end{aligned} \right. \end{equation*} $\end{document} </tex-math></disp-formula></p>
<p>where $ \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\neq 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.</p></abstract>