  The present paper studies the existence of weak solutions for
\begin{equation*} (\mathcal{P}) \left\{\begin{aligned} (-\Delta)^{s_1}_{p_1} u
&=\la f_1\,(x,u,v) +g_1(x,u) \,\mbox{ in }\, \Om, \\ (-\Delta)^{s_2}_{p_2} v
&=\la f_2\,(x,u,v) +g_2(x,v) \,\mbox{ in }\, \Om, \\ u=v &= 0 \,\mbox{in }\,
\Rn \setminus \Om, \\ \end{aligned} \right. \end{equation*} where $\Om \subset
\Rn$ is a smooth bounded domain with smooth boundary, $s_1,s_2 \in (0,1)$,
$1<p_i<\frac{N}{s_i}$, $i=1,2$, $f_i$ and $g_i$ has certain growth assumptions
for $i=1,2$. We prove existence of at least three non negative solutions of
$(\mathcal P)$ under restrictive range of $\lambda$ using variational methods.
As a consequence, we also conclude that a similar result can be obtained when
we consider a more general non local operator $\mathcal L_{\phi_i}$ instead of
$(-\Delta)^{s_i}_{p_i}$ in $(\mathcal P)$.
