We are concerned with the following Lane-Emden system with Neumann boundary conditions: \(\begin{aligned} {\left\{ \begin{array}{ll} -\Delta u_1=|u_2|^{p_\epsilon -1}u_2,\ & in\ \Omega ,\\ -\Delta u_2=|u_1|^{q_\epsilon -1}u_1, \ & in\ \Omega ,\\ \partial _\nu u_1=\partial _\nu u_2=0,\ & on\ \partial \Omega \end{array}\right. } \end{aligned}\) where \(\Omega =B_1(0)\) is the unit ball in \(\mathbb {R}^n\) ( \(n\ge 4\) ) centered at the origin, \(p_\epsilon =p+\alpha \epsilon , q_\epsilon =q+\beta \epsilon \) with \(\alpha ,\beta >0\) and \(\frac{1}{p+1}+\frac{1}{q+1}=\frac{n-2}{n}\) . We construct solutions based on the Lyapunov-Schmidt reduction argument incorporating the zero-average condition by certain symmetries. It’s important to emphasize that we are simultaneously examining two scenarios: when p is greater than \(\frac{n}{n-2}\) and when p is less than \(\frac{n}{n-2}\) . The coupling mechanisms within the Lane-Emden system differ significantly between these two cases, resulting in notable variations in solution behaviors. This distinction also leads to varying research challenges. Presently, there is few literature that takes both of these ranges into consideration when addressing solution construction. Consequently, this aspect constitutes the main feature and new ingredient of our work.
