[proofplan]
We use Darboux's criterion for [Riemann integrability](/page/Riemann%20Integrability). The finite partition witnessing piecewise continuity isolates the only possible discontinuities at finitely many breakpoints. On each open piece, the continuous extension is uniformly continuous, so the oscillation of $f$ on sufficiently small subintervals away from the breakpoints is small. We cover the finitely many breakpoints by intervals of small total length, use boundedness to control their total contribution, and make the Darboux upper sum minus lower sum arbitrarily small.
[/proofplan]
[step:Choose the finite partition and obtain a global bound]
Let
\begin{align*}
a = x_0 < x_1 < \cdots < x_n = b
\end{align*}
be a partition as in the statement. For each $k \in \{1,\dots,n\}$, let
\begin{align*}
g_k: [x_{k-1},x_k] \to \mathbb{R}
\end{align*}
be the continuous extension of $f|_{(x_{k-1},x_k)}$.
Since each $g_k$ is continuous on the compact interval $[x_{k-1},x_k]$, it is bounded there. Define
\begin{align*}
A_k := \sup\{|g_k(x)| : x \in [x_{k-1},x_k]\}
\end{align*}
for each $k \in \{1,\dots,n\}$. Also define the finite breakpoint bound
\begin{align*}
B := \max\{|f(x_j)| : j \in \{0,\dots,n\}\}.
\end{align*}
Now set
\begin{align*}
M := \max\{B,A_1,\dots,A_n,1\}.
\end{align*}
Then $M \geq 1$ and
\begin{align*}
|f(x)| \leq M
\end{align*}
for every $x \in [a,b]$. Indeed, if $x$ is one of the partition points $x_j$, this follows from the definition of $B$; otherwise $x \in (x_{k-1},x_k)$ for some $k$, and $f(x)=g_k(x)$.
[/step]
[step:Separate the breakpoint contribution from the continuous pieces]
Let $\varepsilon > 0$ be given. Define the finite breakpoint set
\begin{align*}
F := \{x_0,x_1,\dots,x_n\}.
\end{align*}
Choose a number $\eta > 0$ such that
\begin{align*}
8M(n+1)\eta < \frac{\varepsilon}{2}.
\end{align*}
For each $j \in \{0,\dots,n\}$, define the interval
\begin{align*}
I_j := [a,b] \cap (x_j-\eta,x_j+\eta).
\end{align*}
Let
\begin{align*}
U := \bigcup_{j=0}^{n} I_j
\end{align*}
be the union of these breakpoint neighborhoods, and define the remaining set
\begin{align*}
K := [a,b] \setminus U.
\end{align*}
The total length of the covering intervals $I_j$ is at most $2(n+1)\eta$ in the following concrete sense: each $I_j$ is contained in an interval of length $2\eta$, so the sum of the lengths of these finitely many covering intervals is at most $2(n+1)\eta$, regardless of overlap and regardless of truncation at $a$ or $b$.
We will classify as bad every partition subinterval that intersects $U$, rather than only those contained in $U$. If a partition subinterval $J$ has length less than $\eta$ and intersects some $I_j$, then every point of $J$ lies in
\begin{align*}
[a,b] \cap (x_j-2\eta,x_j+2\eta).
\end{align*}
Thus the union of all bad subintervals is covered by the finitely many intervals $[a,b] \cap (x_j-2\eta,x_j+2\eta)$, whose total covering length is at most $4(n+1)\eta$. On each bad subinterval the oscillation of $f$ is at most $2M$, so the total possible bad contribution to an upper-minus-lower Darboux sum is bounded by
\begin{align*}
2M \cdot 4(n+1)\eta < \frac{\varepsilon}{2}.
\end{align*}
[/step]
[step:Make the oscillation small away from the breakpoints]
For each $k \in \{1,\dots,n\}$, the function $g_k: [x_{k-1},x_k] \to \mathbb{R}$ is continuous on a compact interval, hence uniformly continuous. Therefore there exists $\delta_k > 0$ such that whenever $u,v \in [x_{k-1},x_k]$ and $|u-v| < \delta_k$, one has
\begin{align*}
|g_k(u)-g_k(v)| < \frac{\varepsilon}{2(b-a)}.
\end{align*}
Define
\begin{align*}
\delta := \min\{\delta_1,\dots,\delta_n,\eta\}.
\end{align*}
Construct a partition $P$ as follows. Start with the finite set consisting of $a$, $b$, every point $x_j$, and every point of $[a,b]$ of the form $x_j-\eta$ or $x_j+\eta$. Discard the points that do not lie in $[a,b]$, discard repetitions, and list the remaining points in increasing order. Refine each resulting interval, if necessary, by inserting finitely many equally spaced points so that every final subinterval has length less than $\delta$. This gives a partition
\begin{align*}
a = y_0 < y_1 < \cdots < y_m = b
\end{align*}
whose mesh satisfies
\begin{align*}
|P| := \max_{1 \leq r \leq m}(y_r-y_{r-1}) < \delta.
\end{align*}
For each subinterval $J_r := [y_{r-1},y_r]$ contained in $K$, the interval $J_r$ lies inside one closed interval $[x_{k-1},x_k]$ and contains no breakpoint. Since $f=g_k$ on the interior of that piece and $J_r$ avoids the breakpoint neighborhoods, the oscillation of $f$ on $J_r$ equals the oscillation of $g_k$ on $J_r$. Because the length of $J_r$ is less than $\delta \leq \delta_k$, we obtain
\begin{align*}
\sup_{x \in J_r} f(x) - \inf_{x \in J_r} f(x) \leq \frac{\varepsilon}{2(b-a)}.
\end{align*}
[guided]
The purpose of removing the set $U$ is to isolate the only points where $f$ may fail to agree with a [continuous function](/page/Continuous%20Function) on a closed interval. Away from $U$, every sufficiently small subinterval sits entirely inside one continuity piece, so [uniform continuity](/page/Uniform%20Continuity) controls the oscillation.
Fix $k \in \{1,\dots,n\}$. The map
\begin{align*}
g_k: [x_{k-1},x_k] \to \mathbb{R}
\end{align*}
is continuous on a compact interval. Hence it is uniformly continuous: for the positive number
\begin{align*}
\frac{\varepsilon}{2(b-a)}
\end{align*}
there exists $\delta_k > 0$ such that for all $u,v \in [x_{k-1},x_k]$ with $|u-v| < \delta_k$,
\begin{align*}
|g_k(u)-g_k(v)| < \frac{\varepsilon}{2(b-a)}.
\end{align*}
We do this for every one of the finitely many pieces and set
\begin{align*}
\delta := \min\{\delta_1,\dots,\delta_n,\eta\}.
\end{align*}
The minimum is positive because it is the minimum of finitely many positive numbers.
Now construct a partition by starting with the finite set consisting of $a$, $b$, every breakpoint $x_j$, and every point of $[a,b]$ of the form $x_j-\eta$ or $x_j+\eta$. After discarding points outside $[a,b]$ and removing repetitions, list the remaining points in increasing order. If one of the resulting intervals has length at least $\delta$, subdivide it into finitely many equally spaced pieces, each of length less than $\delta$. This produces a partition
\begin{align*}
a = y_0 < y_1 < \cdots < y_m = b
\end{align*}
whose mesh satisfies
\begin{align*}
|P| := \max_{1 \leq r \leq m}(y_r-y_{r-1}) < \delta.
\end{align*}
Consider a subinterval
\begin{align*}
J_r := [y_{r-1},y_r]
\end{align*}
that is contained in $K=[a,b]\setminus U$. Because $J_r$ avoids the breakpoint neighborhoods and the partition includes all breakpoints, $J_r$ lies in exactly one continuity piece $[x_{k-1},x_k]$. On the interior of this piece, $f$ agrees with $g_k$, and no endpoint ambiguity from a breakpoint occurs because $J_r$ is away from $F$. Thus the oscillation of $f$ on $J_r$ is controlled by the oscillation of $g_k$ on the same interval.
For any $u,v \in J_r$, we have
\begin{align*}
|u-v| \leq y_r-y_{r-1} < \delta \leq \delta_k.
\end{align*}
Uniform continuity gives
\begin{align*}
|f(u)-f(v)| = |g_k(u)-g_k(v)| < \frac{\varepsilon}{2(b-a)}.
\end{align*}
Taking the supremum over $u \in J_r$ and the infimum over $v \in J_r$ yields
\begin{align*}
\sup_{x \in J_r} f(x) - \inf_{x \in J_r} f(x) \leq \frac{\varepsilon}{2(b-a)}.
\end{align*}
This is the exact estimate needed for the good part of Darboux's criterion: outside the finitely many breakpoint neighborhoods, every subinterval contributes only a small oscillation times its length. The remaining subintervals will be handled as bad intervals. We do not need them to be contained in $U$; it is enough that they intersect $U$, because the mesh bound $|P|<\eta$ places every such interval inside one of the enlarged neighborhoods $[a,b]\cap(x_j-2\eta,x_j+2\eta)$, whose total covering length was bounded earlier.
[/guided]
[/step]
[step:Estimate the Darboux gap for a refined partition]
For each subinterval $J_r=[y_{r-1},y_r]$ of $P$, define
\begin{align*}
M_r := \sup\{f(x):x \in J_r\}
\end{align*}
and
\begin{align*}
m_r := \inf\{f(x):x \in J_r\}.
\end{align*}
The upper and lower Darboux sums of $f$ with respect to $P$ are
\begin{align*}
U(f,P) := \sum_{r=1}^{m} M_r (y_r-y_{r-1})
\end{align*}
and
\begin{align*}
L(f,P) := \sum_{r=1}^{m} m_r (y_r-y_{r-1}).
\end{align*}
Their difference is
\begin{align*}
U(f,P)-L(f,P) = \sum_{r=1}^{m} (M_r-m_r)(y_r-y_{r-1}).
\end{align*}
Split the sum into bad subintervals, those with $J_r \cap U \neq \varnothing$, and good subintervals, those with $J_r \cap U = \varnothing$. Every good subinterval satisfies $J_r \subset K$ by the definition $K=[a,b]\setminus U$, so the oscillation estimate from the previous step applies to the good subintervals.
For bad subintervals, the oscillation bound $M_r-m_r \leq 2M$ holds because $|f|\leq M$ on $[a,b]$. Since $|P|<\delta\leq\eta$, every bad subinterval that intersects $I_j$ is contained in $[a,b]\cap(x_j-2\eta,x_j+2\eta)$. Therefore the union of all bad subintervals is covered by the finitely many intervals $[a,b]\cap(x_j-2\eta,x_j+2\eta)$, whose total covering length is at most $4(n+1)\eta$. Hence the bad contribution is bounded by
\begin{align*}
\sum_{J_r \cap U \neq \varnothing} (M_r-m_r)(y_r-y_{r-1}) \leq 2M \cdot 4(n+1)\eta < \frac{\varepsilon}{2}.
\end{align*}
For good subintervals, the oscillation estimate from the previous step gives
\begin{align*}
\sum_{J_r \cap U = \varnothing} (M_r-m_r)(y_r-y_{r-1}) \leq \frac{\varepsilon}{2(b-a)} \sum_{J_r \cap U = \varnothing}(y_r-y_{r-1}).
\end{align*}
Since the good subintervals lie inside $[a,b]$, their total length is at most $b-a$, and therefore
\begin{align*}
\sum_{J_r \cap U = \varnothing} (M_r-m_r)(y_r-y_{r-1}) \leq \frac{\varepsilon}{2}.
\end{align*}
Combining the two estimates gives
\begin{align*}
U(f,P)-L(f,P) < \varepsilon.
\end{align*}
[/step]
[step:Conclude Riemann integrability by Darboux's criterion]
We have shown that for every $\varepsilon > 0$ there exists a partition $P$ of $[a,b]$ such that
\begin{align*}
U(f,P)-L(f,P) < \varepsilon.
\end{align*}
By Darboux's criterion for Riemann integrability, $f$ is Riemann integrable on $[a,b]$. This proves the theorem.
[/step]
