Proving $S_1-S_2$ is non-positive I have system number (1) as:
$S_1’=-bS_1I_1+(1-c)aI_1$
$I_1’=bS_1I_1-aI_1$
And on the other hand I have system number (2) as:
$S_2’=-bS_2I_2+aI_2$
$I_2’=bS_2I_2-aI_2$
All parameters $a,b,c>0$.
Both systems have the same initial conditions as in $S_1(0)=S_2(0)=s_0>0$ and $I_1(0)=I_2(0)=i_0>0$.
I want to prove that $S_1-S_2$ is non-positive.
What I did so far:

*

*I found $S_2$ explicitly though not sure it is useful as $S_1$ can’t be found explicitly unfortunately.

*Using the initial conditions then I have $S_1=  \frac{(1-c)a}{b}+(s_0- \frac{(1-c)a}{b}$$ )e^{-\int_{0}^{t} bI_1(t) \ dt}$.

*Similarly, I obtained $S_2=  \frac{a}{b}+(s_0- \frac{a}{b}$$ )e^{-\int_{0}^{t} bI_2(t) \ dt}$.

*So from 2 and 3 it is deduced that $(S_1-S_2)(0)=0$ and $$\lim_{t\rightarrow \infty}(S_1-S_2)(t)=-ca/b.$$

*By manipulating the equation for $S_1$ then we have the relationship $e^{-\int_{0}^{t} bI_1(t) \ dt}= \frac{|bS_1(t)-(1-c)a|}{|bs_0-(1-c)a|} $
Despite all that I can’t seem to prove that $(S_1-S_2)(t)$ is negative for $t>0$. I checked this via simulations and indeed it is the case but I’m trying by contradiction and can’t find something that can lead to a contradiction in anything. Help is appreciated!
 A: Contrary to the widespread belief, your system (which looks like an "endemic" version of SIR with acquired immunity, so I naturally wonder whether that is what it really models) can be solved explicitly in the sense that the 3D trajectory $(t,S,I)$ has a parametric representation in terms of two elementary functions and one integral of an elementary function, which allows one to analyze it if not in full, then, at least, to the extent sufficient for answering questions like yours. I'll set $b=1$ (just by scaling time) for convenience.
Replacing $S$ by $S-a$, we see that we need to consider the family of equations
$$
\dot S_\sigma=-I_\sigma(S_\sigma+\sigma),\qquad \dot I_\sigma=S_\sigma I_\sigma
$$
with initial data $S_\sigma(0)=s\in\mathbb R$, $I_\sigma(0)=i>0$ (note that after subtracting $a$ we no longer can assume that $s>0$) and to show that $S_\sigma\le S_0$ for all times if $\sigma>0$.
Clearly, $I_\sigma>0$ for all times and, thereby, $S_\sigma$ goes towards $-\sigma$ in a monotone fashion as $t\to+\infty$ (however, it doesn't need to reach it at $+\infty$, so your analysis of the limit there is flawed as written). So, if $-\sigma\le s\le 0$, $S_0$ goes up, $S_\sigma$ goes down, and the inequality is trivial. Thus we can assume WLOG that $s$ and $s+\sigma$ have the same sign.
Now, given the solution $S_\sigma(t), I_\sigma(t)$,  solve the equation
$$
\frac{dt_{\sigma}}{d\xi}=I_\sigma(t_\sigma(\xi))^{-1}
$$
with the initial condition $t_\sigma(0)=0$.
Note that $S_\sigma+I_\sigma$ is non-increasing for $\sigma>0$ and $S_\sigma$ stays bounded between $s$ and $-\sigma$, so $I_\sigma$ stays bounded as well. Thus the values of $t_\sigma$ cover $[0,+\infty)$ (though $+\infty$ can be reached with finite $\xi$) and, therefore, $t_\sigma$ is a decent time change. Switching to that new variable (i.e., considering $t_\sigma, S_\sigma\circ t_\sigma, I_\sigma\circ t_\sigma$, we get the system
$$
\frac{d}{d\xi} S_\sigma=-(S_\sigma+\sigma)\qquad \frac{d}{d\xi}I_\sigma=S_\sigma,
\qquad \frac{d}{d\xi}t_\sigma=I_\sigma^{-1}\,,
$$
which, taking into account the initial conditions,  readily solves as
$$
S_\sigma(\xi)=-\sigma+(s+\sigma)e^{-\xi},
\\
I_\sigma(\xi)=i-\sigma\xi+(s+\sigma)(1-e^{-\xi}),
\\
t_\sigma(\xi)=\int_0^\xi\frac{d\eta}{i-\sigma\eta+(s+\sigma)(1-e^{-\eta})}\,,
$$
where the meaningful values of $\xi$ are those for which $I_\sigma>0$ on $[0,\xi]$.
Due to the simultaneous monotonicity of $S_0$ and $S_\sigma$, the inequality we want to prove is equivalent to the statement that if $S_\sigma(\xi_\sigma)=S_0(\xi_0)$, then $t_\sigma(\xi_\sigma)\le t_0(\xi_0)$ if $s,s+\sigma>0$ and $t_\sigma(\xi_\sigma)\le t_0(\xi_0)$ if $s,s+\sigma<0$.
To prove this statement, let us define $\lambda=\lambda(\eta)\ge 0$ by
$$
S_\sigma(\lambda)=-\sigma+(s+\sigma)e^{-\lambda}=se^{-\eta}=S_0(\eta)
$$
as long as $S_0(\eta)$ makes sense and the equation is solvable.
We derive from this definition
$$
e^{\lambda-\eta}=-\frac{\sigma}se^{\lambda}+\frac{s+\sigma}{s}
\\
d\lambda=\frac{s}{s+\sigma}e^{\lambda-\eta}\,d\eta=\left[-\frac{\sigma}{s+\sigma}e^{\lambda}+1\right]\,d\eta
\\
(s+\sigma)(1-e^{-\lambda})=s(1-e^{-\eta})\,.
$$
Using all these identities together with $\xi_\sigma=\lambda(\xi_0)$ and making the change of variable $\eta\mapsto\lambda(\eta)$ in the integral defining $t_\sigma(\xi_\sigma)$, we see that we have to compare
$$
t_\sigma(\xi_\sigma)=\int_0^{\xi_0}\frac
{\left[-\frac{\sigma}{s+\sigma}e^{\lambda}+1\right]\,d\eta}
{i-\sigma\lambda+s(1-e^{-\eta})}
$$
and
$$
t_0(\xi_0)=\int_0^{\xi_0}\frac{d\eta}{i+s(1-e^{-\eta})}\,.
$$
We shall just compare the integrands. Multiplying by the denominators (which should stay non-negative throughout $[0,\xi_0]$ for the assumption $S_\sigma(t_\sigma)=S_0(t_0)$ to be possible) and cancelling equal terms we get to compare
$$
-\frac{\sigma}{s+\sigma}e^{\lambda}(i+s(1-e^{-\eta})\text{ and }-\sigma\lambda
$$
on $[0,\xi_0]$.
Case 1: $s,s+\sigma<0$ (increasing). Then $LHS\ge 0\ge RHS$.
Case 2: $s,s+\sigma>0$ (decreasing). Then (since $i>0$),
$$
LHS\le -\frac{\sigma}{s+\sigma}e^{\lambda}s(1-e^{-\eta})=
-\sigma e^{\lambda}(1-e^{-\lambda})
\\
=-\sigma(e^\lambda-1)\le -\sigma\lambda=RHS
$$
and we are done.
