Hint (borrowing from my answer to the homogeneous problem): eliminating $z$ and $y\,$, respectively, between the linear equations:
$$
\begin{cases}
(a_1c_2-a_2c_1)x+(b_1c_2-b_2c_1)y = d_1c_2−d_2c_1 \\
(a_1b_2-a_2b_1)x+(b_2c_1-b_1c_2)z = d_1b_2−d_2b_1
\end{cases}
\\ \iff
\quad
\begin{cases}
(a_1c_2-a_2c_1)x - (d_1c_2−d_2c_1) = -(b_1c_2-b_2c_1)y \\
(a_1b_2-a_2b_1)x - (d_1b_2−d_2b_1) = -(b_2c_1-b_1c_2)z
\end{cases}
$$
Multiplying the latter equations:
$$
\begin{align}
\big((a_1c_2-a_2c_1)x - (d_1c_2−d_2c_1)\big)\big((a_1b_2-a_2b_1)x - (d_1b_2−d_2b_1)\big) &= (b_1c_2-b_2c_1)(b_2c_1-b_1c_2)yz \\ &= (b_1c_2-b_2c_1)(b_2c_1-b_1c_2)C
\end{align}
$$
The above is a quadratic in $x$ (in the general case, unless some of the coefficients are $\,0$), and solvability in reals depends on the sign of the discriminant thereof.
The method of characteristics looks like the right way to solve this. Along paths that satisfy ${\rm d}x_i/{\rm d}t = f_i(\vec{x})$, one finds $u(\vec{x}(t))$ evolves according to ${\rm d}u/{\rm d}t = -c u$. If the path terminates at $\partial\Omega$, then $u(x) = 0$ along the whole path. This leads to our first necessary condition for the existence of a nonzero solution:
(1) $\exists$ path $\vec{x}(t)$ satisfying ${\rm d}x_i/{\rm d}t = f_i(\vec{x})$ with origin and terminus (limits as $t \rightarrow \pm\infty$) in the interior of $\Omega$.
For a continuous $u(\vec{x})$, the value of $u(\vec{x}(t))$ cannot diverge when $t \rightarrow \pm\infty$. Excepting a set of measure zero, all paths $\vec{x}(t)$ start at a repulsor and end at an attractor (rather than, say, a saddle point). Two more necessary conditions for the existence of a nonzero solution are therefore:
(2) $c < 0$ at $\vec{x}(-\infty)$
(3) $c > 0$ at $\vec{x}(+\infty)$
Except for a set of measure zero, we can probably assume these inequalities are strict, i.e. $c < 0$ and $c > 0$, respectively (convergence is possible for $c = 0$ but not guaranteed, depending on derivative terms). With the strict inequalities, conditions (1-3) are also sufficient for nonzero solutions $u(\vec{x})$ to exist. That can be seen as follows:
Starting with a point $\vec{x}_0$ along the path $\vec{x}(t)$, define a size-$\epsilon$ cross section (orthogonal to the streamlines of ${\rm d}x_i/{\rm d}t = f_i(\vec{x})$) and posit that $u(\vec{x})$ varies smoothly from $u(x_0) = 1$ to $u = 0$ at the boundaries of the cross section. The value of $u(\vec{x})$ along the "past" and "future" of this cross section is obtained by propagating along the characteristics using ${\rm d}u/{\rm d}t = -c u$. All these characteristics originate from the same repulsor (where $u = 0$) and terminate at the same attractor (also where $u = 0$). Fill in the rest of $\Omega$ with the null solution $u = 0$. Thus we have constructed a nonzero, continuous-valued solution to the PDE.
There are a bunch of singular edge cases where the necessary and sufficient conditions don't coincide, i.e. if $\lVert f \rVert = u = 0$ at the same point (fixable by rescaling $f$ and $u$), if $\lVert f\rVert = 0$ over an open subset of $\Omega$, if $\lVert f\rVert = 0$ on the boundary $\partial\Omega$, if $c = 0$ at $\vec{x}(\pm\infty)$. In the space of possible functions $(\vec{f}, u)$, these singular cases only occur in a set of measure zero, so are not very interesting. Almost everywhere, conditions (1-3) are both necessary and sufficient.
Putting this another way, we can say (almost everywhere) that the zero solution is unique if:
$\forall$ paths $\vec{x}(t)$ satisfying ${\rm d}x_i/{\rm d}t = f_i(\vec{x})$ with origin and terminus in the interior of $\Omega$,
$c > 0$ at $\vec{x}(-\infty)$ or $c < 0$ at $\vec{x}(+\infty)$.
Coming back to your condition $c^* < 0$: Note that $\partial_i f^i < 0$ at attractors (this always holds, regardless of whether it's a node, limit cycle, toroid, chaotic attractor, etc.). Therefore, if $c^* < 0$ on $\Omega$, it follows that $c = c^* + \partial_i f^i < 0$ at all of the attractors. Therefore, the second condition above is always satisfied when $c^* < 0$. The condition above is the more general sufficient (and necessary) condition for uniqueness (with the caveats noted above).
Since any dynamical system can be represented by ${\rm d}x_i/{\rm d}t = f_i(\vec{x})$ and dynamical systems can be really, really complicated, the general condition can be hard to work with, so more specific conditions like $c^* < 0$ might be more useful.
Also, defining the value of $c$ is tricky when the attractor / repulsor isn't a point. Taking the average over limit cycles is straightforward, chaotic attractors less so (ergodic theory).
Best Answer
Differentiate the first equation w.r.t. $x_2$ and second equation w.r.t. $x_1$. The equality of mixed partials holds for smooth solutions, so that $f_{x_2} = g_{x_1}$ must be satisfied. Now, integrate the first differential equation w.r.t. $x_1$ as follows: $$ X(x_1,x_2) - X(a,x_2) = \int_{a}^{x_1} f(\xi,x_2)\, d\xi \, . $$ Then, injecting this expression in the second differential equation gives $$ \int_{a}^{x_1} f_{x_2}(\xi,x_2)\, d\xi + X_{x_2}(a,x_2) = g(x_1,x_2)-g(a,x_2) + X_{x_2}(a,x_2) = g(x_1,x_2) \, , $$ which is naturally satisfied. Of course, if we integrate the second equation first, we get $$ X(x_1,x_2) - X(x_1,b) = \int_{b}^{x_2} g(x_1,\eta)\, d\eta \, . $$ Both expressions are linked through the condition $f_{x_2} = g_{x_1}$. You may also have a look at this post and related links.
For a boundary condition of the form $X(0,x_2)=0$, we find $$ X(x_1,x_2) = \int_0^{x_1} f(\xi,x_2)\, d\xi $$ by setting $a=0$.