My attempt:
I want to take advantage of the fact that I can "section off" infinitely many of the discontinuities that are close to zero. I'll be using the following technique in my proof: if for every $\epsilon > 0$, there is a partition $P_\epsilon$ such that $U(f,P_\epsilon)-L(f,P_\epsilon)<\epsilon$, then $f$ is Riemann integrable. I'll divide my solution into steps:
- The function is bounded, say by $B \geq |f(x)|$. Choose $\alpha$ such that $0 < \alpha < \frac{\epsilon}{6B}$. Let $P_1$ simply be the partition $\{0,\alpha\}$. Then, using the notation $M_k = \sup_{x\in [x_{k-1},x_k]} f(x)$ and $m_k = \inf_{x\in [x_{k-1},x_k]} f(x)$, we have
$$ U(f,P_1)-L(f,P_1) = \sum_{k}(M_k-m_k)(x_k-x_{k-1}) = (M_1-m_1)(\alpha-0) \leq 2B \frac{\epsilon}{6B} = \frac{\epsilon}{3}$$
- Now consider the interval $[\alpha,1]$. There are only finitely many discontinuities in here; say $N$ of them. We'll only consider part of this interval. We can break it up into disjoint closed intervals which contain the remaining discontinuous points. Let $\beta < \frac{\epsilon}{6NB}$. Let $P_2$ be the partition $\{x_0, ..., x_{2N}\}$ such that the intervals $[x_i, x_{i+1}]$ are of length $\beta$ and contain the discontinuities as their midpoint. I.e, they are of the form $[1/n-\beta/2, 1/n+\beta/2]$. If any of these intervals overlap, just choose $\beta$ to be even smaller so that they do not. Now,
$$U(f,P_2)-L(f,P_2) = \sum_{k=1}^{2N}(M_k-m_k)(x_k-x_{k-1}) \leq 2B (N\beta) < 2B\frac{\epsilon}{6NB}N = \frac{\epsilon}{3}$$
Finally, consider the partition $\{y_0, ..., y_P\}$ consisting of the "remaining" points; i.e. those such that $$ [0,\alpha] \cup \left( \bigcup_{k=1}^{2N} [x_{k-1},x_k]\right) \cup \left( \bigcup_{k=1}^{P} [y_{k-1},y_k]\right) = [0,1]$$ $f$ is continuous and bounded on these intervals, so it is Riemann integrable on these intervals. Thus we can break them into subintervals $P_3 := \{y_0, ..., y_L\}$, (with $L \geq P$), such that $$U(f,P_3)-L(f,P_3) = \sum_{k}(M_k-m_k)(x_k-x_{k-1}) < \frac{\epsilon}{3}$$
Now we can combine the first 3 steps: Let $P_\epsilon = P_1 \cup P_2 \cup P_3$, so that $P_\epsilon$ is a partition for $[0,1]$ satisfying $$U(f,P_\epsilon)-L(f,P_\epsilon) < \frac{\epsilon}{3} + \frac{\epsilon}{3} + \frac{\epsilon}{3} = \epsilon $$
Does this proof look correct?
The intuition is right, and the formalization is almost right, but there are just some small technical details which could be improved on (but atleast you've clearly conveyed what you're trying to do, so these small details are easily identified, and with "a little" more effort, they can be fixed).
You should choose something like $0<\alpha < \min\left(\frac{1}{2}, \frac{\epsilon}{6B}\right)$, just to ensure that $0<\alpha<1$ (which is what you've implicitly been assuming).
I realize that for $(2)$, your idea is to sandwich each discontinuity point $\frac{1}{n}$ inside some interval $[x_{i-1}, x_i]$ of very small length, however, the estimates you wrote down are false. Note that it is the subsequent intervals which have small length, so for example $[x_0,x_1], [x_2, x_{3}], [x_4,x_5], [x_6,x_7]$ etc have small length, but the intervals in between like $[x_1,x_2], [x_5,x_6], [x_7,x_8]$ need not have small length (see the picture below).
Because of this, the estimate \begin{align} \sum_{k=1}^{2N}(M_k-m_k)(x_k - x_{k-1}) \leq (2B) (N\beta) \end{align} is false. You should probably have terms like $(x_{2k+1} - x_k)$ appearing in the sum. Also, another small detail is that you say you want to create the partition $P_2$ such that the discontinuities are at the midpoint of the intervals of length $\beta$.
Well, this isn't always possible, because suppose $f$ is discontinuous at the endpoint $1$. How are you supposed to keep $1$ at the midpoint of the interval when $[0,1]$ is the domain of definition?
There may be other small details I may have missed, but the point is this: you're trying to do too many things too explicitly in a single proof, each of which is linked back to a previous step. So, if you provide an incorrect proof for one of the steps, the rest of the proof is in jeopardy. I suggest you proceed as follows:
Now, once you establish this lemma, you're good to go: let $\epsilon > 0$, and $B>0$ a bound on $f$. Then, choose $0<\alpha< \min\left( 1, \frac{\epsilon}{4B}\right)$. Now, since $f$ has finitely many discontinuities on $[\alpha,1]$ by the lemma above, it follows $f$ is integrable on $[\alpha,1]$. So, there a partition $P'$ of $[\alpha,1]$ such that \begin{align} U(f,P') - L(f,P') < \frac{\epsilon}{2}. \end{align} Now, by taking $P= \{0\}\cup P'$ it follows (since $\alpha\in P'$ by definition) a simple calculation shows that \begin{align} U(f,P)-L(f,P) < (2B)(\alpha) + \frac{\epsilon}{2} < \epsilon. \end{align}
Note that by treating the case of a single discontinuity first as a lemma, there are very few details to keep track of, so you're less likely to make a mistake in the proof (even minor ones), and it's easier to keep track of all the cases. Next, the induction step is pretty trivial in this situation. If you structure your proof in this manner, there is less opportunity for mistakes here.
Anyway, if you want to continue working on your original approach, then what you essentially have to do is transcribe the above picture properly into words (which is definitely doable, but I find it too cumbersome to get all the $\epsilon,\alpha,\beta$ and the various partitions right).
Just an FYI, there's a much stronger result: Lebesgue's Criterion for Riemann-integrability.