In this paper, we consider the system − Δ u = λ ( v + 1 ) p , − Δ v = γ ( u + 1 ) θ -\Delta u =\lambda (v+1)^p,\;\;-\Delta v = \gamma (u+1)^\theta on a smooth bounded domain Ω \Omega in R N \mathbb {R}^N with Dirichlet boundary condition u = v = 0 u=v=0 on ∂ Ω . \partial \Omega . Here λ , γ \lambda ,\gamma are positive parameters and 1 > p ≤ θ 1 > p \le \theta . Let x 0 x_0 be the largest root of the polynomial H ( x ) = x 4 − a m p ; 16 p θ ( p + 1 ) ( θ + 1 ) ( p θ − 1 ) 2 x 2 + 16 p θ ( p + 1 ) ( θ + 1 ) ( p + θ + 2 ) ( p θ − 1 ) 3 x a m p ; − 16 p θ ( p + 1 ) 2 ( θ + 1 ) 2 ( p θ − 1 ) 4 . \begin{align*} H(x) = x^4 - &\frac {16p\theta (p+1)(\theta +1)}{(p\theta -1)^2}x^2 + \frac {16p\theta (p+1)(\theta +1)(p+\theta +2)}{(p\theta -1)^3}x\\ &-\frac {16p\theta (p+1)^2(\theta +1)^2}{(p\theta -1)^4}. \end{align*} We show that the extremal solutions associated to the above system are bounded provided N > 2 + 2 x 0 . N>2+2x_0. This improves the previous work by Craig Cowan (2015). We also prove that if N ≥ 2 + 2 x 0 , N\geq 2+2x_0, then the singular set of any extremal solution has Hausdorff dimension less than or equal to N − ( 2 + 2 x 0 ) . N-(2+2x_0).