The bounds change due to $u$. The original lower bound was $y=0$, so
$$u(\text{lower bound})=4y-3y^2+6y^3+1=4(0)-3(0)^2+6(0)^3+1=1$$
which is the lower bound. Do the same process and you will get the upper bound.
As @TheCount mentions, this is because when one uses the substitution, the end result is a function of $u$. You then substitute $u=f(y)$ back into it, and then apply the bounds.
Or you can apply the bounds first and not worry about re-substitution later.
First, it is a sufficient condition.
As long as $f$ is bounded on $[0,1]$, the upper and lower sums corresponding to arbitrary partitions are bounded. Let $\mathcal{P}$ denote the set of all partitions of $[0,1].$ Consequently, the sets $\{L(P,f): P \in \mathcal{P}\}$ and $\{U(P,f): P \in \mathcal{P}\}$ are bounded, and this guarantees the existence of
$$\underline{\int}_0^1 f(x) \, dx = \sup_{P \in \mathcal{P}}\, L(P,f), \\ \overline{\int}_0^1 f(x) \, dx = \inf_{P \in \mathcal{P}}\, U(P,f) , $$
which are called the lower and upper integrals.
Given any regular partition $D_n$ we have
$$L(D_n,f) \leqslant \underline{\int}_0^1 f(x) \, dx \leqslant \overline{\int}_0^1 f(x) \, dx \leqslant U(D_n,f).$$
The central inequality follows because for any partitions $P$ and $Q$ we have $L(P,f) \leqslant U(Q,f)$ (take a common refinement of the partitions to show this) and, thus $\sup_{P \in \mathcal{P}} \,L(P,f) \leqslant \inf_{Q \in \mathcal{P}} \,U(Q,f)$.
Hence,
$$0 \leqslant \overline{\int}_0^1 f(x) \, dx - \underline{\int}_0^1 f(x) \, dx \leqslant U(D_n,f) - L(D_n,f).$$
The right-hand side converges to $0$ as $n \to \infty$, by hypothesis, which implies that $f$ is integrable since we must have
$$\underline{\int}_0^1 f(x) \, dx = \overline{\int}_0^1 f(x) \, dx, $$
where the common value of lower and upper integrals is by definition the value of the integral.
To show it is a necessary condition, consider
$$\left|U(D_n,f) - L(D_n,f) \right| \leqslant \left|U(D_n,f) - \int_0^1 f(x) \, dx \right| + \left|L(D_n,f) - \int_0^1 f(x) \, dx \right|.$$
The two terms on the RHS go to zero as $n \to \infty$. This is a consequence of the equivalent condition for integrability where for arbitrary Riemann sums corresponding to tagged partitions we have
$$\tag{*}\int_0^1 f(x) \, dx = \lim_{\|P\| \to 0} S(P,f).$$
Here $\|P\| = \max_{1 \leqslant j \leqslant n} (x_j - x_{j-1})$ is the norm of the partition $P = (x_0,x_1, \ldots, x_n)$ and, clearly, $\|D_n\| \to 0$ if and only if $n \to \infty$.
It takes a bit of effort to prove the equivalence of $(*)$ to the definition of the Riemann integral in terms of partition refinement or the Darboux approach. It has been shown a number of times on this site including here.
Best Answer
The process is that the function in the lower limit needs to be less than the function in the upper limit over the interval of integration. In your case, $y^2 \leq y+2$ for all $y$ in $[-1,2]$.
As for how to determine that inequality, you could do it by algebra, but the quickest way is probably with a graph. You know that $x=y^2$ is a parabola with its vertex at the origin, and which opens to the right. The line $y=x-2$ passes through $(0,-2)$ and $(2,0)$.
A few more comments:
You may think that you can just try it with the functions in one order, and if you get a negative number, just drop the minus sign and make it positive. This may work, but only if the graphs don't cross each other in the interval over which you're integrating.
In multivariable calculus, you really need graphing skills. Trying to do it all by just pattern matching and algebra doesn't help you gain any intuition.
You asked which “equation” should be in the lower and upper limits of integration, but you mean which function. For instance $y=x+2$ is an equation, but when you set up the integral it had to be manipulated into $x=y-2$, because you needed a function of $y$.