If $f$ is integrable on $[a,b]$ and $[b,c]$ then $f$ is integrable on $[a,c]$ and $$\int_a^c f(x) dx = \int_a^b f(x) dx + \int_b^c f(x) dx .$$
Let's define three partitions $P,Q,R$ on $[a,c]$, $[a,b]$, $[b,c]$ respectively: $$P = \{x_0, x_1, ..., x_k, ..., x_n\}$$ $$Q = \{x_0, x_1, ..., x_k\}$$ $$R = \{x_k, ..., x_n\}$$ where $x_0 = a, x_k = b, x_n= c$
As the function $f$ is integrable on both $[a,b]$ and $[b,c]$, for any $\epsilon/2>0$ there is a $\delta>0$ such that $$\left| S(Q,f) - \int_a^b f(x) dx\right|<\epsilon/2 $$ and $$\left| S(R,f) - \int_b^c f(x) dx \right|<\epsilon/2 $$ Adding both inequalities $$\left| S(Q,f) + S(R,f) - \int_a^b f(x) dx - \int_b^c f(x) dx \right|< \epsilon$$ Should we assume that there no uniform subdivision of the partitions, given $x_j\leq d_j \leq x_{j-1}$ we have $$ \sum_{j=1}^{k} f(d_j)(x_j-x_{j-1}) + \sum_{j=k}^{n} f(d_j) (x_j-x_{j-1}) \geq \sum_{j=1}^{n} f(d_j) (x_j - x_{j-1})$$ $$ S(R,f) +S(Q,f) \geq S(P,f)$$
taking the limits where $||P|| \rightarrow 0$ or $n \rightarrow \infty$, (given $H'$ a more refined partition of $H$, we have $S(H')\leq S(H)$) $$\lim\limits_{n \rightarrow \infty} S(R,f) + \lim\limits_{n \rightarrow \infty} S(Q,f) \leq \lim\limits_{n \rightarrow \infty} S(P,f)$$ $$ \int_a^b f(x) dx + \int_b^c f(x) dx \leq \int_c^a f(x) dx$$
It shows that $$\left| S(P,f) - \int_a^c f(x) dx \right| \leq \left| S(Q,f) + S(R,f) - \int_a^b f(x) dx - \int_b^c f(x) dx \right|< \epsilon$$
It follows that for any $\epsilon>0$ there is a $\delta>0$ such that $$\left| S(P,f) - \int_a^c f(x) dx \right| < \epsilon, ||P||<\delta$$
Hence, $f$ is by definition integrable on $[a,c]$
Here is the proof I did. Is my argumentation valid, or does it need some improvement? Beside, I have difficulty in understanding/applying the theorem $S(P)-s(P)<\epsilon$, how can it be used here? Much appreciated for your help/input.
You should be careful to work consistently within one of the equivalent frameworks for defining the Riemann integral (Darboux sums or Riemann sums) and where convergence to the integral is based on partition mesh or refinement.
You started along one path correctly and then fell short of the mark. Let's proceed as you did starting with the definition of the integral over $[a,c]$, that for every $\epsilon > 0$ there is a $\delta > 0$ such that if $P$ is a partition with $\|P\| < \delta$ then
$$\left|S(P,f;[a,c]) - \int_a^cf(x)\, dx\right| < \epsilon.$$
I have introduced the notation $S(P,f;[a,c])$ to indicate the interval on which the partition and sum are defined.
Given $\epsilon > 0$ we have $\delta_1, \delta_2$ such that $\|P_i\| < \delta_i$ for $i = 1,2$ implies
$$\left|S(P_1,f;[a,b]) - \int_a^bf(x)\, dx\right| < \epsilon/3, \\ \left|S(P_2,f;[b,c]) - \int_b^cf(x)\, dx\right| < \epsilon/3, $$
Let $M$ be an upper bound for $|f|$ and take $\delta = \min(\delta_1, \delta_2, \epsilon/6M)$.
Let $P = (x_0,x_1, \ldots, x_{j-1},x_j, \ldots , x_n)$ be a partition of $[a,c]$ with $\|P\| < \delta.$
In general $P$ will not include the point $b$ (if it does the proof is easy). I will assume WLOG that $x_{j-1} < b < x_j$.
In this case,
$$S(P,f;[a,c]) = \sum_{k=1}^{j-1}f(\xi_k)(x_k - x_{k-1}) + f(\xi_j)(x_j - x_{j-1}) + \sum_{k=j}^{n}f(\xi_k)(x_k - x_{k-1}),$$
where $\xi_k \in [x_{k-1},x_k]$. Again WLOG I assume $x_{j-1} \leqslant \xi_j < b < x_j$.
Now we can write $f(\xi_j)(x_j - x_{j-1}) = f(\xi_j)(b- x_{j-1}) + (f(\xi_j) - f(b))(x_j - b) + f(b)(x_j - b), $ and
$$\begin{align}S(P,f;[a,c]) &= \underbrace{\sum_{k=1}^{j-1}f(\xi_k)(x_k - x_{k-1}) + f(\xi_j)(b- x_{j-1})}_{S(P_1,f;[a,b])}\\ &+ (f(\xi_j) - f(b))(x_j - b) \\ &+ \underbrace{\sum_{k=j}^{n}f(\xi_k)(x_k - x_{k-1}) + f(b)(x_j- b)}_{S(P_2,f;[b,c])}\end{align}$$
Hence, since $|f(\xi_j) - f(b)|(x_{j} - b) < 2M\delta < \epsilon/3$,
$$\left|S(P,f;[a,c]) - \left(\int_a^b f(x) \, dx + \int_b^c f(x) \, dx\right) \right| \\ \leqslant \left|S(P_1,f;[a,b]) - \int_a^b f(x) \, dx \right|+ |f(\xi_j) - f(b)|(x_{j} - b) + \left|S(P_2,f;[b,c]) - \int_b^c f(x) \, dx \right| \\ < \epsilon. $$
Thus, $f$ is integrable over $[a,c]$ and
$$\int_a^c f(x) \, dx = \int_a^b f(x) \, dx + \int_b^c f(x) \, dx. $$