How to prove that continuous functions are Riemann-integrable?

First note that a precise formulation of your question is:

How do you prove that every continuous function on a closed bounded interval is Riemann (not Darboux) integrable?

You can find a proof in Chapter 8 of these notes.

Here is a rough outline of this handout:

I. I introduce the ("definite") integral axiomatically. One of the axioms is that the set of integrable functions on $[a,b]$ should contain all the continuous functions.
II. I prove that the Fundamental Theorem of Calculus follows (easily) from the axioms.
III. I introduce Riemann integrable functions (which are exactly what you wrote above) and verify that the class of Riemann integrable functions on $[a,b]$ satisfies the axioms of I. In particular:
IV: I prove that every continuous function is Riemann integrable.

Later I talk about the Darboux integral and how it compares to the Riemann integral. But it was an intentional decision to present the Riemann integral first. This is what students are expecting from their previous courses, and it is not so bad to work with, at least for a while.

Riemann-Darboux integral:

For the Riemann-Darboux integral it is easier than with your definition (which is equivalent in $\mathbf R$).

The function is uniformly continuous on $[a,b]$ (why?). This means that if $\epsilon > 0$ is given we can find $\delta = \delta(\epsilon) > 0$ such that $|x - y| < \delta$ implies $|f(x) - f(y)| < \frac{\epsilon}{2(b - a)}$. Now let $P_\epsilon$ be a partition with norm $\|P_\epsilon\| < \delta$. Now for $P$ finer than $P_\epsilon$ we have

$$M_k(f) - m_k(f) \leq \frac{\epsilon}{2(b - a)}.$$

Where $M_k$ is the supremum in $[x_{k - 1}, x_k]$ and $m_k$ is the infimum. Multiply this inequality with $\Delta x_k$ and sum to get

$$U(P, f) - L(P, f) \leq \frac{\epsilon}{2(b - a)} \sum_{k = 1}^n \Delta x_k = \frac\epsilon2 < \epsilon.$$

Fine, so what is left is to prove the same thing for your definition of the integral (not so easy) or proving the equivalence between the two (not so hard).

Riemann integral:

There is a problem with the above approach if we are in a general Banach space (why?), so we must resort to the normal Riemann integral.

Let $f:[a,b] \to E$ be continuous where $E$ is a Banach space. Given $\epsilon > 0$ let $\delta$ be such that $$\text{if } |x - y| < \delta \text{ then } |f(x) - f(y)| < \frac{\epsilon}{b - a}.$$

Let $P$ and $P'$ be partitions of $[a,b]$ with norm smaller than $\delta$. Let $c$ and $c'$ be the choices of points in each interval of $P$ and $P'$ respectively. We want to estimate $|S(P, c) - S(P', c')|$ where $S(P, c)$ is the Riemann sum associated with $P$ and $c$. WLOG let $P \subset P'$. (If $P = P'$ then $$|S(P, c) - S(P, c')| \leq \sum |f(c_i) - f(c_i')| \Delta x_i \leq \epsilon.$$

Now suppose that $P'$ is obtained from $P$ by inserting one point (split one interval) for example say we insert $x_j'$ with $x_j \leq x_j' \leq x_{j + 1}$. In this case the partition size will not increase. WLOG assume that for $i \neq j$ we have $x_i' = x_i$ and that $c_j = x_j'$ and that $x_j'$ is also selected as the points in the intervals $[x_j, x_j']$ and $[x_j', x_{j + 1}]$. Now $S(P, c) - S(P', c') = 0$. We can repeat this process for another refinement.

This will give us the result (why?).

Appendix to Chapter 13 in: CALCULUS, by M. Spivak.
Yes, uniform continuity is the key.