Introduce binaries $\delta_1$and $\delta_2$ satisfying $\delta_1+\delta_2=1$. Now you have two cases, $\delta_1=1$, $c_1^Tx=0$, $x_1+x_2-x_3\geq 0$ and $\delta_2=1$, $c_2^Tx=0$, $x_1+x_2-x_3\leq -\epsilon$. Standard big-M modelling and you have something like
$$ -M(1-\delta_1) \leq c_1^Tx \leq M(1-\delta_1)\\
-M(1-\delta_2) \leq c_2^Tx \leq M(1-\delta_2)\\
x_1+x_2-x_3 \geq -M(1-\delta_1)\\
x_1+x_2-x_3 \leq -\epsilon + M(1-\delta_2)
Of course, in practice you use different big-M constants on each constraint and try to make them as small as possible.
$\delta_1+\delta_2 =1$ is redundant as the two sets are disjoint, but it never hurts to add simple constraints on binaries.