Inhomogeneous Wave Equation on Bounded Domain

As you said, $S_1$ is the linear propagator of the equation. That is to say that $S_1(t) \begin{bmatrix}\phi\\\psi\end{bmatrix} $ is by definition the solution of the homogeneous equation at time $t$ with initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix} $. So, you can write it explicitly once you solve the homogeneous equation for all initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix}$, and you do that through the fourier-series method. So, the right hand side of $$S_1(t) \begin{bmatrix}\phi\\\psi\end{bmatrix} = \sum_{n=1}^\infty \bigl[A_n\cos(nt) + B_n\sin(nt)\bigr]\sin(nx)$$ is just what you obtain, and the constants $A_n,\, B_n,\, C_n,\, D_n$ are explicit too, as they are essentially the fourier coefficients of the initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix} $ (you can find them in the end of page 29).

Note that $S_1(t) \begin{bmatrix}\phi\\\psi\end{bmatrix} $ is not a vector. $S_1(t)$ does not act like a square matrix or something like that. $S_1(t) \begin{bmatrix}\phi\\\psi\end{bmatrix} $ is just the solution of the equation with initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix} $, so it’s a scalar. Note also that the first formula holds for an arbitrary couple of initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix} $ and for an arbitrary time $t\geq 0$.

Then, the last equation is a consequence of the first one, it’s just that your initial data is $\begin{bmatrix}0\\f(s)\end{bmatrix} $, so you have to change the coefficients accordingly (you have to use the fourier coefficients of $f(s)$) and the time you take is not $t$, but $t-s$. You just have to use the same formula we talked about before, as it describes the way the operator $S_1(t)$ acts on a couple of initial data $\begin{bmatrix}\phi\\\psi\end{bmatrix} $, and you have to specialize that formula to the initial data $\begin{bmatrix}0\\f(s)\end{bmatrix} $.

I hope this is clear now. Sorry for the repetitions