  We study positive solutions to the steady state reaction diffusion systems of
the form: \begin{equation} \left\{\begin{array}{ll} -\Delta u = \lambda
f(v)+\mu h(u), & \Omega,\\ -\Delta v = \lambda g(u)+\mu q(v),& \Omega,\\
\frac{\partial u}{\partial \eta}+\sqrt[]{\lambda +\mu}\, u=0,&
\partial\Omega,\\ \frac{\partial v}{\partial \eta}+\sqrt[]{\lambda +\mu}\, v=0,
& \partial\Omega,\\ \end{array}\right. \end{equation} where ${\lambda,\mu>0}$
are positive parameters, ${\Omega}$ is a bounded in $\mathbb{R}^{N}$$(N>1)$
with smooth boundary ${\partial \Omega}$, or ${\Omega=(0,1)}$, ${
\frac{\partial z}{\partial \eta} }$ is the outward normal derivative of $z$.
Here $f, g, h, q\in C^{2} [0,r)\cap C[0,\infty)$ for some $r>0$. Further, we
assume that $f, g, h$ and $q$ are increasing functions such that $f(0) = g(0) =
h(0) = {q}(0) = 0$, $f^\prime(0), g^\prime(0), h^\prime(0), q^\prime(0) > 0$,
and $\lim\limits_{s\to \infty}\frac{f(M g(s))}{s}=0$ for all $M>0$. Under
certain additional assumptions on $f, g, h$ and $ q$ we prove our existence and
multiplicity results. Our existence and multiplicity results are proved using
sub-super solution methods.
