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 spectral 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.