Here is an alternative approach to solving the original equation.

The desired equation subject to the constraint $U \ge 0$ (elementwise nonnegativity) can be written as

\begin{equation*}
U = C + A^TUA,\qquad U \ge 0.
\end{equation*}

Define now the iteration:

\begin{equation*}
U_{k+1} = [C + A^TU_kA]_+,\qquad k=0,1,\ldots
\end{equation*}
where $[\cdot]_+ \equiv \max\lbrace 0, \cdot\rbrace$ denotes **projection** onto the nonnegative orthant. Assuming the slightly stronger condition that the operator norm $\|A\| < 1$, we can show that the above iteration converges to the desired solution.

$\newcommand{\Gc}{\mathcal{G}}$
**Claim.** Let $\Gc$ denote the map $U\mapsto [C + A^TUA]_+$. Then, $\Gc$ is a strict contraction in the spectral norm.

*Proof.* Let $U, V$ be two matrices. Then,
\begin{eqnarray*}
\|\Gc(U)-\Gc(V)\| &=& \|[C + A^TUA]_+ - [C + A^TVA]_+\|\\
&\le& \|C + A^TUA - C - A^TVA]\|\\
&\le& \|A^T(U-V)A\|\\
&\le& \|A^T\| \|A\| \|U-V\| \\
&\le& \gamma\|U-V\|,\qquad \text{where}\ \gamma=\|A\|^2 < 1.
\end{eqnarray*}
The first inequality follows as projection is a *non-expansive* operator, the rest follow from obvious properties of the spectral norm.

Thus, $\Gc$ is a strict contraction. So by the Banach contraction theorem the iteration $U_{k+1}=\Gc(U_k)$ has a unique solution. This solution is nonnegative by construction and satisfies the desired equation.