  We concern a family $\{(u_{\varepsilon},v_{\varepsilon})\}_{\varepsilon > 0}$
of solutions of the Lane-Emden system on a smooth bounded convex domain
$\Omega$ in $\mathbb{R}^N$ \[\begin{cases} -\Delta u_{\varepsilon} =
v_{\varepsilon}^p &\text{in } \Omega,\\ -\Delta v_{\varepsilon} =
u_{\varepsilon}^{q_{\varepsilon}} &\text{in } \Omega,\\ u_{\varepsilon},\,
v_{\varepsilon} > 0 &\text{in } \Omega,\\ u_{\varepsilon} = v_{\varepsilon} =0
&\text{on } \partial\Omega \end{cases}\] for $N \ge 4$,
$\max\{1,\frac{3}{N-2}\} < p < q_{\varepsilon}$ and small \[\varepsilon :=
\frac{N}{p+1} + \frac{N}{q_{\varepsilon}+1} - (N-2) > 0.\] This system appears
as the extremal equation of the Sobolev embedding $W^{2,(p+1)/p}(\Omega)
\hookrightarrow L^{q_{\varepsilon}+1}(\Omega)$, and is also closely related to
the Calder\'on-Zygmund estimate. Under the a natural energy condition
\[\sup_{\varepsilon > 0} \left(\|u_{\varepsilon}\|_{W^{2,{p+1 \over
p}}(\Omega)} + \|v_{\varepsilon}\|_{W^{2,{q_{\varepsilon}+1 \over
q_{\varepsilon}}}(\Omega)}\right) < \infty,\] we prove that the multiple
bubbling phenomena may arise for the family
$\{(u_{\varepsilon},v_{\varepsilon})\}_{\varepsilon > 0}$, and establish a
detailed qualitative and quantitative description. If $p < \frac{N}{N-2}$, the
nonlinear structure of the system makes the interaction between bubbles so
strong, so the determination process of the blow-up rates and locations is
completely different from that of the classical Lane-Emden equation. If $p \ge
\frac{N}{N-2}$, the blow-up scenario is relatively close to (but not the same
as) that of the classical Lane-Emden equation, and only one-bubble solutions
can exist. Even in the latter case, the standard approach does not work well,
which forces us to devise a new method. Using our analysis, we also deduce a
general existence theorem valid on any smooth bounded domains.
