We examine regularity of the extremal solution of nonlinear nonlocal eigenvalue problem Lu = λF (u, v) in Ω, Lv = γG (u, v) in Ω,with an integro-differential operator, including the fractional Laplacian, of the formwhen J is a nonnegative measurable even jump kernel. In particular, we consider jump kernels of the form of J(y) = a(y/|y|) |y| n+2s where s ∈ (0, 1) and a is any nonnegative even measurable function in L 1 (S n−1 ) that satisfies ellipticity assumptions. We first establish stability inequalities for minimal solutions of the above system for a general nonlinearity and a general kernel. Then, we prove regularity of the extremal solution in dimensions n < 10s and n < 2s + 4s p∓1 [p + p(p ∓ 1)] for the Gelfand and Lane-Emden systems when p > 1 (with positive and negative exponents), respectively. When s → 1, these dimensions are optimal. However, for the case of s ∈ (0, 1) getting the optimal dimension remains as an open problem. Moreover, for general nonlinearities, we consider gradient systems and we establish regularity of the extremal solution in dimensions n < 4s. As far as we know, this is the first regularity result on the extremal solution of nonlocal system of equations.