Geometry provides insight and classical inequalities afford easy access to rigor.
Geometric solution
We know, from the geometry of least squares, that $\mathbf{\bar{x}} = (\bar{x}, \bar{x}, \ldots, \bar{x})$ is the orthogonal projection of the vector of data $\mathbf{x}=(x_1, x_2, \ldots, x_n)$ onto the linear subspace generated by the constant vector $(1,1,\ldots,1)$ and that $\sigma_x$ is directly proportional to the (Euclidean) distance between $\mathbf{x}$ and $\mathbf{\bar{x}}.$ The non-negativity constraints are linear and distance is a convex function, whence the extremes of distance must be attained at the edges of the cone determined by the constraints. This cone is the positive orthant in $\mathbb{R}^n$ and its edges are the coordinate axes, whence it immediately follows that all but one of the $x_i$ must be zero at the maximum distances. For such a set of data, a direct (simple) calculation shows $\sigma_x/\bar{x}=\sqrt{n}.$
Solution exploiting classical inequalities
$\sigma_x/\bar{x}$ is optimized simultaneously with any monotonic transformation thereof. In light of this, let's maximize
$$\frac{x_1^2+x_2^2+\ldots+x_n^2}{(x_1+x_2+\ldots+x_n)^2} = \frac{1}{n}\left(\frac{n-1}{n}\left(\frac{\sigma_x}{\bar{x}}\right)^2+1\right) = f\left(\frac{\sigma_x}{\bar{x}}\right).$$
(The formula for $f$ may look mysterious until you realize it just records the steps one would take in algebraically manipulating $\sigma_x/\bar{x}$ to get it into a simple looking form, which is the left hand side.)
An easy way begins with Holder's Inequality,
$$x_1^2+x_2^2+\ldots+x_n^2 \le \left(x_1+x_2+\ldots+x_n\right)\max(\{x_i\}).$$
(This needs no special proof in this simple context: merely replace one factor of each term $x_i^2 = x_i \times x_i$ by the maximum component $\max(\{x_i\})$: obviously the sum of squares will not decrease. Factoring out the common term $\max(\{x_i\})$ yields the right hand side of the inequality.)
Because the $x_i$ are not all $0$ (that would leave $\sigma_x/\bar{x}$ undefined), division by the square of their sum is valid and gives the equivalent inequality
$$\frac{x_1^2+x_2^2+\ldots+x_n^2}{(x_1+x_2+\ldots+x_n)^2} \le \frac{\max(\{x_i\})}{x_1+x_2+\ldots+x_n}.$$
Because the denominator cannot be less than the numerator (which itself is just one of the terms in the denominator), the right hand side is dominated by the value $1$, which is achieved only when all but one of the $x_i$ equal $0$. Whence
$$\frac{\sigma_x}{\bar{x}} \le f^{-1}\left(1\right) = \sqrt{\left(1 \times (n - 1)\right)\frac{n}{n-1}}=\sqrt{n}.$$
Alternative approach
Because the $x_i$ are nonnegative and cannot sum to $0$, the values $p(i) = x_i/(x_1+x_2+\ldots+x_n)$ determine a probability distribution $F$ on $\{1,2,\ldots,n\}$. Writing $s$ for the sum of the $x_i$, we recognize
$$\eqalign{
\frac{x_1^2+x_2^2+\ldots+x_n^2}{(x_1+x_2+\ldots+x_n)^2} &= \frac{x_1^2+x_2^2+\ldots+x_n^2}{s^2} \\
&= \left(\frac{x_1}{s}\right)\left(\frac{x_1}{s}\right)+\left(\frac{x_2}{s}\right)\left(\frac{x_2}{s}\right) + \ldots + \left(\frac{x_n}{s}\right)\left(\frac{x_n}{s}\right)\\
&= p_1 p_1 + p_2 p_2 + \ldots + p_n p_n\\
&= \mathbb{E}_F[p].
}$$
The axiomatic fact that no probability can exceed $1$ implies this expectation cannot exceed $1$, either, but it's easy to make it equal to $1$ by setting all but one of the $p_i$ equal to $0$ and therefore exactly one of the $x_i$ is nonzero. Compute the coefficient of variation as in the last line of the geometric solution above.
