We examine regularity of the extremal solution of nonlinear nonlocal eigenvalue problem {Lu = lambda F(u, v) in Omega, Lv = gamma G(u, v) in Omega, u, v = 0 on R-n \ Omega, with an integro-differential operator, including the fractional Laplacian, of the form L(u(x)) = lim(epsilon -> 0) integral(Rn\B epsilon(x)) [u(x) - u(z)]J(z - x)dz, when J is a nonnegative measurable even jump kernel. In particular, we con- sider jump kernels of the form of J(y) = a(y/vertical bar y vertical bar)/vertical bar y vertical bar(n+2s) where s is an element of (0, 1) and a is any nonnegative even measurable function in L-1(Sn-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 +root 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 is an element of (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.