  We study the behavior near the origin of $C^2$ positive solutions $u(x)$ and
$v(x)$ of the system
  $0\le -\Delta u \le (\frac{1}{|x|^\alpha}* v)^\lambda$
  $0\le -\Delta v \le (\frac{1}{|x|^\beta}* u)^\sigma$
  in $B_2(0)\setminus\{0\} \subset R^n$, $n\ge 3$, where $\lambda,\sigma \ge 0$
and $\alpha,\beta\in (0,n)$.
  A by-product of our methods used to study these solutions will be results on
the behavior near the origin of $L^1(B_1(0))$ solutions $f$ and $g$ of the
system
  $0 \le f(x) \le C(|x|^{2-\alpha} + \int_{|y|<1}\frac{ g(y)
dy}{|x-y|^{\alpha-2}} )^\lambda$
  $0 \le g(x) \le C(|x|^{2-\beta} + \int_{|y|<1}\frac{ f(y)
dy}{|x-y|^{\beta-2}} )^\sigma$
  for $0<|x|<1$ where $\lambda,\sigma \ge 0$ and $\alpha, \beta\in (2,n+2)$.
