Proof of particular solution to a DE in the form of $y''+a_1y'+a_2y = A\cos (bx) + B\sin (bx)$ In the book "Differential Equations and Linear Algebra" by Stephen W. Goode (2nd edition), he presents a theorem 2.4.1 (pg 128) where in part 1 he states that the DE
$$y''+a_1y'+a_2y = A\cos (bx) + B\sin (bx)$$
where $a_1$, $a_2$, $A$, $B$, and $b$ are constants, has a particular solution of the form
$$y_p(x)=A_0\cos(bx)+B_0\sin(bx)$$
if $r=ib$ is not a root of the auxiliary equation. Goode leaves the proof of this theorem as an exercise to the reader so I started with the simple substitution of the $y_p$ into the DE as such
$$-b^2(A_0\cos(bx)+B_0\sin(bx))+a_1b(-A_0\sin(bx)+B_0\cos(bx))+a_2(A_0\cos(bx)+B_0\sin(bx))=A\cos(bx)+B\sin(bx)$$
but I cannot figure out how to manipulate this equation to prove that $y_p$ is a solution if $P(-ib)\not=0$
 A: Let $r_1,r_2$ be the solutions of $r^2+a_1r+a_2=0\iff r^2-(r_1+r_2)r+r_1r_2=0$
You can rewrite your equation as a system of two first order equations:
$y''-(r_1+r_2)y'+r_1r_2y=\overbrace{(y''-r_1y')}^{u'}-r_2\overbrace{(y'-r_1y)}^u=f(x)$
$\begin{cases}y'-r_1y=u\\u'-r_2u=f(x)\end{cases}$
So the homegeneous solution of $u$ is $u_h=Ce^{r_2x}$ that you can also rewrite in sincos form when $r_2$ is a complex root.
From there it all depends whether $b$ is equal to $r_2$ or not. When applying the variation of the constant (i.e. $C$) method to find a particular solution. It'll be $A_0\cos(bx)+B_0\sin(bx)$ in the case $b\neq r_2$, else you will have linear coefficients, not constant ones.
The same applies then when you solve the second equation $y'-r_1y=u=u_h+u_p$, this time we need $b\neq r_1$ and $r_1\neq r_2$ to keep the particular solution in the same form.

Edit:
Let try to find a particular solution to RHS $f(x)=A\cos(bx)+B\sin(bx)$.


*

*method 1: variation of the constant

$\require{cancel}u=Ce^{rx}\implies u'-ru=\cancel{Cre^{rx}}+C'e^{rx}-\cancel{Cre^{rx}}=f(x)=\frac{A-iB}2e^{ibx}+\frac{A+iB}2e^{-ibx}$
So you get $$C'=\frac{A-iB}2e^{(ib-r)x}+\frac{A+iB}2e^{-(ib+r)x}$$
You can see your discussion about $r\neq\pm ib$ appearing because we converted back $\sin$ and $\cos$ to their respective exponential expression.
From there in the case $r\neq\pm ib$ we can integrate $e^{(ib-r)x}$ to $\frac{e^{(ib-r)x}}{ib-r}$ and same for the other term. When you remultiply by $e^{rx}$ to get back $u_p=Ce^{rx}$ you are left only with constant coefficients and $e^{\pm ibx}$ which in turn can be converted to $\cos(bx)$ and $\sin(bx)$.
Though the coefficients in question being complex, the calculation is really tedious because you need to rationalize the denominator, then convert back to sincos, this is a lot of work.
So instead we consider this step as the theoretical justification that we will end up with $A_0\cos(bx)+B_0\sin(bx)$ and we prefer to use direct identification below.


*

*method 2: direct identification

The theoretical justification above lead to this more general formulation


*

*your homogeneous equation has root $b$ with multiplicity $m$ (convention $m=0$ if $b$ is not a root)

*the full equation has a RHS of the form $P(x)e^{bx}$ with $P$ polynomial.

 Then you need to search for a particular solution in the form
$Q(x)e^{bx}$ with $Q$ polynomial and $$\deg(Q)=\deg(P)+m$$
Although since the homogeneous solution will already have vanishing
terms $(C_0+C_1x+\cdots+C_{m-1}x^{m-1})e^{bx}$, you can ignore them in
the polynomial Q.
In case of a linear combination of such terms in the RHS, you can also
search for a linear combination of particular solutions for each.
Note: in the special case of $RHS = P(x)$, consider the root $b=0$ since $e^{0x}=1$, and the same rule applies.

Therefore since we need to find only $1$ particular solution, we can just apply this recipe and go for identification:
$u_p=A_0\cos(bx)+B_0\sin(bx)$
$u_p'-ru_p=(B_0b-A_0r)\cos(bx)-(A_0b+B_0r)\sin(bx)$
And you get to solve $\begin{cases}B_0b-A_0r = A\\A_0b-B_0r = -B\end{cases}$
The final result is as below:
$$\overbrace{-\frac{Ar+Bb}{r^2+b^2}}^{A_0}\cos(bx)+\overbrace{\frac{Ab-Br}{r^2+b^2}}^{B_0}\sin(bx)$$


We have done it with the system of ODE of first order to illustrate, but the general formulation above works of course directly for the higher order linear ODE with constant coefficients.
So having $y''+a_1y'+a_2y=A\cos(bx)+B\sin(bx)$
You set the characteristic equation $r^2+a_1r+a_2=0$ find roots $r_1,r_2$, compare them to $\pm ib$ (the exponential form of sincos) and determine the $m$ above.
In our case $r\neq\pm ib\implies m=0$ and we get to search for a particular solution of the same degree (i.e. constant coefficients) than the initial RHS.
i.e. $RHS=A\cos(bx)+B\sin(bx)\implies u_p=A_0\cos(bx)+B_0\sin(bx)$
