Solution (source code)

= Solution

Let $E$ be the underlying <Banach space>. By the <Laplace-transform formula for a semigroup resolvent>,
$$
R_\lambda f=\int_0^\infty e^{-\lambda s}P_sf\,ds,\qquad
\|R_\lambda\|\leq\int_0^\infty e^{-\lambda s}\,ds=\frac1\lambda.
$$
These are <Bochner integrals>, defined for every $f\in E$ and $\lambda>0$. The <contraction semigroup> bound and <semigroup property> give
$$
R_\lambda R_\mu f
=\int_0^\infty\int_0^\infty e^{-\lambda s-\mu t}P_{s+t}f\,dt\,ds.
$$
Interchanging and regrouping these <Bochner integrals> is justified by the integrable bound $e^{-\lambda s-\mu t}\|f\|$, whose integral is $\|f\|/(\lambda\mu)$. For $\lambda\ne\mu$, set $v=s+t$ and evaluate the inner scalar integral:
$$
\begin{aligned}
R_\lambda R_\mu f
&=\int_0^\infty\left(\int_0^ve^{-\lambda s-\mu(v-s)}\,ds\right)P_vf\,dv\\
&=\frac1{\mu-\lambda}\int_0^\infty(e^{-\lambda v}-e^{-\mu v})P_vf\,dv\\
&=\frac{R_\lambda f-R_\mu f}{\mu-\lambda}.
\end{aligned}
$$
Since $f$ was arbitrary, the <resolvent identity> follows:
$$
\boxed{R_\lambda-R_\mu=(\mu-\lambda)R_\lambda R_\mu.}
$$
For $\lambda=\mu$, both sides vanish, so that case holds as well.