We study systems of equations of the form $X_1 = f_1(X_1, \ldots, X_n), \ldots, X_n = f_n(X_1, \ldots, X_n)$ where each $f_i$ is a polynomial with nonnegative coefficients that add up to~$1$. The least nonnegative solution, say~$\mu$, of such equation systems is central to problems from various areas, like physics, biology, computational linguistics and probabilistic program verification. We give a simple and strongly polynomial algorithm to decide whether $\mu=(1,\ldots,1)$ holds. Furthermore, we present an algorithm that computes reliable sequences of lower and upper bounds on~$\mu$, converging linearly to~$\mu$. Our algorithm has these features despite using inexact arithmetic for efficiency. We report on experiments that show the performance of our algorithms.