If
\(x=x_0\text{,}\) the formula is immediate. If
\(x \lt x_0\text{,}\) then the same argument applied on the interval with endpoints
\(x\) and
\(x_0\) gives the result after reversing the limits of integration. So it is enough to treat the case
\(x_0 \lt x\text{.}\)
For each \(k=0,1,\dots,n\text{,}\) define
\begin{equation*}
R_k(x;x_0):=\frac{1}{k!}\int_{x_0}^x f^{(k+1)}(t)(x-t)^k\,dt.
\end{equation*}
\begin{equation*}
R_0(x;x_0)=\int_{x_0}^x f'(t)\,dt = f(x)-f(x_0).
\end{equation*}
\begin{equation*}
u(t)=f^{(k)}(t)
\qquad \text{and} \qquad
v(t)=\frac{(x-t)^k}{k!}.
\end{equation*}
Then \(u'(t)=f^{(k+1)}(t)\) and
\begin{equation*}
v'(t)=-\frac{(x-t)^{k-1}}{(k-1)!}.
\end{equation*}
Therefore
\begin{equation*}
R_k(x;x_0)
=
\left[\frac{f^{(k)}(t)(x-t)^k}{k!}\right]_{t=x_0}^{t=x}
+
\frac{1}{(k-1)!}\int_{x_0}^x f^{(k)}(t)(x-t)^{k-1}\,dt.
\end{equation*}
Since \((x-x)^k=0\text{,}\) this becomes
\begin{equation*}
R_k(x;x_0)
=
-\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k + R_{k-1}(x;x_0).
\end{equation*}
Starting from \(R_0(x;x_0)=f(x)-f(x_0)\) and applying the recurrence successively for \(k=1,2,\dots,n\text{,}\) we obtain
\begin{equation*}
R_n(x;x_0)
=
f(x)-\sum_{k=0}^n \frac{f^{(k)}(x_0)}{k!}(x-x_0)^k
=
f(x)-P_n(x;x_0).
\end{equation*}
Rearranging gives the desired formula.