Distribution of $\max(u,v)$, where $u$ and $v$ exponential random variables with mean $1$. Could anyone tell me where's the mistake in calculating the distribution of $u$ and $v$, two independent r.v. following an exponential distribution with parameter $1$?
My calculation :
\begin{align*}
P[\max(u,v) \leq c]&=\int_{0}^{+\infty}\int_{0}^{+\infty}1_{\{\max(u,v) \leq c\}}e^{-u}e^{-v}du dv\\
&=\int_{0}^{+\infty} \left[ \int_{0}^{v}1_{\{\max(u,v)\leq c\}}e^{-u}+\int_{v}^{+\infty}1_{\{\max(u,v)\leq c\}}e^{-u}\right]e^{-v}du dv\\
&=\int_{0}^{+\infty}  \int_{0}^{v}1_{\{v\leq c\}}e^{-u}e^{-v}du dv+\int_{0}^{+\infty}\int_{v}^{+\infty}1_{\{u\leq c\}}e^{-u}e^{-v}du dv\\
&=\int_{0}^{c}  \int_{0}^{v}e^{-u} du\, e^{-v} dv+\int_{0}^{+\infty}\int_{v}^{c}e^{-u} du \,e^{-v} dv
\end{align*}
 A: Does it count if I say taking a double integral? If so, read on. If not, disregard.
We have
\begin{align*}
P(\max\{U,V\}\leq c)&=P(U\leq c, V\leq c)\\
&=P(U\leq c)P( V\leq c)\tag 1\\&=
(1-e^{-c})(1-e^{-c})\tag 2\\
&=(1-e^{-c})^2
\end{align*}
where $(1)$ is true by independence, and $(2)$ is the well known cdf of an exponential random variable.
A: The answer by @probablyme is correct (+1): 
If the max of two values is less than $c,$ then they both have to be
less than $c.$ I'm sure of the equations below because they're in a book
I co-authored, and they have been vetted by numerous reviewers and
(by now) students.
If $X = \max(U, V),$ then
its CDF is $F_X(x) = (1 - e^{-x})^2,$ for $x > 0.$
Also, upon differentiation, its PDF is $f_X(x) = 2e^{-x} - 2e^{-2x}.$
The following simulation in R of 10,000 realizations of $X$ illustrates
the PDF by comparison with a histogram of the simulated
distribution (left panel) and the CDF by comparison with its empirical CDF (right panel). In the right panel, the ECDF is
plotted in dots that appear to merge at places into a heavy line
and the green CDF plot runs very directly through those points.
 u = rexp(10^4);  v = rexp(10^4);  x = pmax(u,v)
 par(mfrow=c(1,2))
   hist(x, prob=T, col="wheat")
     curve(2*exp(-x) - 2*exp(-2*x), 0, 12, col="darkgreen", lwd=2, add=T)
   plot.ecdf(x)
     curve((1 - exp(-x))^2, 0, 5, col="green", add=T)
 par(mfrow=c(1,2))


