Proof

We prove both implications:

Assume, \(f\) is Riemann integrable.
Then, by definition, its Riemann upper and lower integrals exist and are equal.
Therefore, \(f\) is bounded by some step functions \(\phi:[a,b]\mapsto\mathbb R\) and \(\psi:[a,b]\mapsto\mathbb R\) with \[\phi\le f\le \psi.\]
From the definition of supremum and infimum it follows for every \(\epsilon > 0\): \[\int_a^b\psi(x)dx-\int_a^b\phi(x)dx\le \epsilon.\]

Now, assume that \(f\) is bounded and that for every \(\epsilon > 0\) there exist some step functions \(\phi:[a,b]\mapsto\mathbb R\) and \(\psi:[a,b]\mapsto\mathbb R\) with \[\phi\le f\le \psi\quad\quad( * )\] and \[\int_a^b\psi(x)dx-\int_a^b\phi(x)dx\le \epsilon.\]
This means that, as \(\epsilon\) tends to \(0\), the Riemann lower integral of the step function \(\phi\) and the Riemann upper integral of the step function \(\psi\) tend to each other: \[\lim_{\epsilon\to 0}\int_{a~*}^{b}\phi(x)dx=\lim_{\epsilon\to 0}\int_{a}^{b~*}\psi(x)dx.\quad\quad ( * * )\]
On the other hand, from \(( * )\), from the monotony of Riemann integrals for step functions, and from the fact that convergence preserves lower and upper bounds, it follows. \[\lim_{\epsilon\to 0}\int_{a~*}^{b}\phi(x)dx\le \int_{a~*}^{b}f(x)dx \le \int_{a}^{b}f(x)dx\le \int_{a}^{b~*}f(x)dx \le \lim_{\epsilon\to 0}\int_{a}^{b~*}\psi(x)dx.\]
Together with \(( * * )\), this means that \(f\) is Riemann integrable, because in fact, \[ \int_{a~*}^{b}f(x)dx = \int_{a}^{b}f(x)dx= \int_{a}^{b~*}f(x)dx.\]
∎

Thank you to the contributors under CC BY-SA 4.0!

Forster Otto: "Analysis 1, Differential- und Integralrechnung einer Veränderlichen", Vieweg Studium, 1983