In this paper we prove existence and regularity of weak solutions for the following system \(\begin{aligned} {\left\{ \begin{array}{ll} & -\text{ div }\Bigg (\bigg (\Vert \nabla u\Vert ^{p}_{L^{p}}+\Vert \nabla v\Vert ^{p}_{L^{p}}\bigg )|\nabla u|^{p-2}\nabla u\Bigg ) + g(x,u,v)=f \ \ \ \text{ in } \ \Omega ; \\ & -\text{ div }\Bigg (\bigg (\Vert \nabla u\Vert ^{p}_{L^{p}}+\Vert \nabla v\Vert ^{p}_{L^{p}}\bigg )|\nabla v|^{p-2}\nabla v\Bigg ) = h(x,u,v) \ \ \ \ \text{ in } \ \Omega ; \\ & u=v=0 \ \text{ on } \ \partial \Omega . \end{array}\right. } \end{aligned}\) where \(\Omega \) is an open bounded subset of \(\mathbb {R}^N\) , \(N>2\) , \(f\in L^m(\Omega )\) , where \(m>1\) and g, h are two Carathéodory functions, which may be non monotone. We prove that under appropriate conditions on g and h, there is gain of Sobolev and Lebesgue regularity for the solutions of this system.