  We consider a system of differential equations with nonlinear Steklov
boundary conditions, related to the fractional problem $$(-\Delta)^s u_i =
f_i(x,u_i) - \beta u_i^p \sum_{j\neq i} a_{ij} u_j^p,$$ where $i = i,\dots, k$,
$s\in(0,1)$, $p>0$, $a_{ij}>0$ and $\beta>0$. When $k=2$ we develop a
quasi-optimal regularity theory in $C^{0,\alpha}$, uniformly w.r.t. $\beta$,
for every $\alpha < \alpha_{\mathrm opt}=min(1,2s)$; moreover we show that the
traces of the limiting profiles as $\beta\to+\infty$ are Lipschitz continuous
and segregated. Such results are extended to the case of $k\geq3$ densities,
with some restrictions on $s$, $p$ and $a_{ij}$.
