Given pdf of $I$ and $R$ (both $I$ and $R$ are independent RV's), how to find cdf of $W =I^2R$?

Where,

$$
\begin{align}
f_I(i)&=6i(1-i), &0 \leq i \leq 1 \\
f_R(r)&=2r, &0 \leq r\leq 1.
\end{align}
$$

The simplest and surest way to compute the distribution density or probability of a random variable is often to compute the means of functions of this random variable. In the case at hand, one wants to write $\mathrm E(g(W))$ as
$$
\color{blue}{\mathrm E(g(W))=\int g(w)f(w)\mathrm{d}w},
$$
for every bounded measurable function $g$, then one can be sure that $f$ is the density of the distribution of $W$. So, in a way, the functions $g$ play the role of a dummy variable and one wants the equality above to hold for every $g$.

Naturally $W=I^2R$ hence $\mathrm E(g(W))$ is a priori a double integral, but one can be sure that a change of variable will save the day. So, applying the definitions,
$\mathrm E(g(W))=\mathrm E(g(I^2R))$ and
$$
\mathrm E(g(I^2R))=\iint g(x^2y)\cdot[0\leqslant x\leqslant 1]\cdot6x(1-x)\cdot[0\leqslant y\leqslant 1]\cdot2y\cdot\mathrm{d}x\mathrm{d}y,
$$
where, for every property $\mathfrak{A}$, Iverson bracket $[\mathfrak{A}]$ denotes $1$ if $\mathfrak{A}$ holds and $0$ otherwise. 

(Begin of rant: no, I do not like to put the limits of the domain of integration on the integral signs, and yes, I prefer to use the notation $[\mathfrak{A}]$ or its cousin $\mathbb{1}_\mathfrak{A}$ because they are more systematic and, at least to me, less error prone. End of rant.) 

Now, what change of variable? For one of the two new variables, we want $w=x^2y$, of course. For the other, a sensible choice (but not the only one) is $z=x$. The new domain is $0\leqslant w\leqslant z^2\leqslant 1$ and the Jacobian is given by $\mathrm{d}x\mathrm{d}y=z^{-2}\mathrm{d}w\mathrm{d}z$, hence
$$
\mathrm E(g(W))=\int g(w)[0\leqslant w\leqslant 1]\left(\int [w\leqslant z^2\leqslant 1]\cdot6z(1-z)(2wz^{-2})z^{-2}\mathrm{d}z\right)\mathrm{d}w.
$$
By identification, the density $f(w)$ is the quantity enclosed by the parenthesis, that is, for every $0\leqslant w\leqslant1$,
$$
f(w)=\int [w\leqslant z^2\leqslant 1]6z(1-z)(2wz^{-2})z^{-2}\mathrm{d}z=12w\int_{\sqrt{w}}^1 z^{-3}(1-z)\mathrm{d}z,
$$
Finally,
$$
\color{red}{f(w)=6(1-\sqrt{w})^2\cdot[0\leqslant w\leqslant1]}.
$$
