The paper is concerned with the Neumann boundary problems: \begin{equation*} \left\{\begin{array}{l} -\Delta _{p} u+\lambda_{1} u^{p-1}=|u|^{p^{*}-2}u+\frac{\alpha}{p^{*}}|u|^{\alpha-2} |v|^{\beta} u,\,\,\,\,\,x\in \Omega\\ -\Delta _{p} v+\lambda_{2} v^{p-1}=|v|^{p^{*}-2}v+\frac{\beta}{p^{*}}|u|^{\alpha} |v|^{\beta -2} v,\,\,\,\,\,x\in \Omega\\ \frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0,\,\,\,\,\,\,\,\,x\in \partial \Omega\\ u>0,v>0\,\,\,\,\,\,\,x\in \Omega. \end{array}\right. \end{equation*} where $\Delta _{p} u= div(|\nabla u|^{p-2}\nabla u)$, $\alpha, \beta >1, \alpha +\beta=p^{*}=\frac{Np}{N-p}$, and $\Omega$ is a bounded domain with a smooth $C^2$ boundary in $\mathbb{R}^N(N\geq 3)$, and $n$ denotes the unit outward normal to $\partial \Omega$. For various parameters $\alpha, \beta, \lambda_{1}, \lambda_{2}, N,p$, we dicuss the existence of positive solutions for the system.
