Here is my proof.
Theorem: If $f$ is Riemann integrable on $[a,b],c\in[a,b]$, and $g(x)=f(x)$ except for finitely many points $c_1,\cdots,c_k$ in $[a,b]$, then $g$ is Riemann integrable on $[a,b]$, and $\int_{a}^{b}f(x)dx=\int_{a}^{b}g(x)dx.$
Suppose another function $g$ such that $g(x)=f(x)$ except for finitely many points $c_1,\cdots,c_k$. Since $f$ is Riemann integrable, for any $\epsilon>0$, there exists $\delta=\frac{\epsilon}{4n\sum_{i=1}^{k}|g(c_i)-f(c_i)|}$ such that for any partition $P=(a=x_0,\cdots,x_n=b)$ with $\left\lVert P\right\rVert<\delta$, we have $\left|R(f,P)-\int_{a}^{b}f\right|<\frac{\epsilon}{2}.$ Now, \begin{align*} \left|R(g,P)-\int_{a}^{b}f\right|&=\left|R(g,P)-R(f,P)+R(f,P)-\int_{a}^{b}f\right|\\ &\leq\left|R(g,P)-R(f,P)\right|+\left|R(f,P)-\int_{a}^{b}f\right|\\ &<\left|R(g,P)-R(f,P)\right|+\frac{\epsilon}{2}\\ &<\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon.\\ \end{align*}where \begin{align*} |R(g,P)-R(f,P)|&=\left|\sum_{i=1}^{n}g(x_i*)(x_i-x_{i-1})-\sum_{i=1}^{n}f(x_i*)(x_i-x_{i-1})\right|\\ &=\left|\sum_{i=1}^{n}(g(x_i^*)-f(x_i^*))(x_i-x_{i-1})\right|\\ &<\sum_{i=1}^{k}|g(c_i)-f(c_i)|(2n\delta)\\ &=2n\sum_{i=1}^{k}|g(c_i)-f(c_i)|\left(\frac{\epsilon}{4n\sum_{i=1}^{k}|g(c_i)-f(c_i)|}\right)\\ &<\frac{\epsilon}{2} \end{align*}
Is my proof of this theorem correct?