Solve $y''+y=p(x)$ where $p$ is a polynomial function The equation is $$y'' + y=p(x)$$ Let $p(x)=\sum\limits_{k=0}^n a_k x^k$
The homogeneous solution to the equation is $y_h = a\cos(x)+b\sin(x)$. I can look for a particular solution of the form  $y_p=\sum\limits_{k=0}^n b_k x^k$. 
We get $$\sum\limits_{k=0}^{n-2} ((k+2)(k+1)b_{k+2}+b_k)x^k+b_{n-1}x^{n-1}+b_nx^n=\sum\limits_{k=0}^n a_kx^k$$  Hence, if $0\le k \le n-2$, $(k+2)(k+1)b_{k+2}+b_k=a_k$ and $b_{n-1}=a_{n-1}$, $b_n=a_n$
I'm trying to show that $y_p=\sum\limits_{n=0}^{\infty}(-1)^n p^{(2n)}(x)$ from here, how do I proceed?
 A: As you stated, you have $b_n=a_n$ and $b_{n-1}=a_{n-1}$. You also have that, for $0 \le k \le n - 2$ that
$$(k+2)(k+1)b_{k+2}+b_k=a_k \implies b_k = a_k - (k+2)(k+1)b_{k+2} \tag{1}\label{eq1A}$$
Thus, you can start from $k = n - 2$ and work down to $k = 0$ to get each value of $b_k$ in turn. For example with the next $3$ values of $b_k$, you have
$$\begin{equation}\begin{aligned}
b_{n-2} & = a_{n-2} - (n)(n-1)b_{n} \\
& = a_{n-2} - (n)(n-1)a_{n}
\end{aligned}\end{equation}\tag{2}\label{eq2A}$$
$$\begin{equation}\begin{aligned}
b_{n-3} & = a_{n-3} - (n-1)(n-2)b_{n-1} \\
& = a_{n-3} - (n-1)(n-2)a_{n - 1}
\end{aligned}\end{equation}\tag{3}\label{eq3A}$$
$$\begin{equation}\begin{aligned}
b_{n-4} & = a_{n-4} - (n-2)(n-3)b_{n-2} \\
& = a_{n-3} - (n-2)(n-2)(a_{n-2} - (n)(n-1)a_{n})
\end{aligned}\end{equation}\tag{4}\label{eq4A}$$
With these coefficients $b_k$ determined, as Mick's question comment suggests, you can then calculate the coefficients of each $x^k$ of what you're trying to show, i.e., $y_p=\sum\limits_{n=0}^{\infty}(-1)^n p^{(2n)}(x)$, to see whether or not they match the determined $b_k$ coefficients.
A: Using the Laplace transform
$$
s^2Y(s) + Y(s) = P(s)+ sy_0+y'_0
$$
and then
$$
Y(s) = \frac{P(s)}{s^2+1}+\frac {sy_0}{s^2+1}+\frac{y'_0}{s^2+1}
$$
and after inversion
$$
y(x) = \sin(x)*p(x) +y_0\cos x+y'_0 \sin x
$$
where $ \sin(x)*p(x) =\int_0^x\sin(x-\zeta)p(\zeta)d\zeta = \int_0^x\sin(\zeta)p(x-\zeta)d\zeta$ is the convolution operation.
NOTE
$$
\sin(x-\zeta)p(\zeta) = \sum_{k=0}^n a_k \zeta^k(\sin x\cos\zeta-\cos x\sin\zeta)
$$ 
the integrals
$$
\int_0^x\zeta^k\cos\zeta d\zeta,\ \ \ \int_0^x\zeta^k\sin\zeta d\zeta
$$
can be handled as a sub-product from
$$
\int_0^x e^{i\zeta}\zeta^k d\zeta
$$
which obeys to the recursion
$$
I_k(x) = -i\left(x^k e^{i x}-k I_{k-1}(x)\right)
$$
for $k > 0$
A: You know that your particular solution is a polynomial. So there is some $m$ so that $y_p^{(2m)}=0$ along with all higher order derivatives. Now apply the binomial formula
$$
A^m-B^m=(A-B)(A^{m-1}+A^{m-2}B+...+B^{m-1})
$$
with $A=D^2$ and $B=-1$ to the differential equation
$$
D^{2m}y_p-(-1)^my_p=(D^{2m-2}-D^{2m-4}+D^{2m-6}-...+(-1)^{m-1})\underbrace{(D^2+1)y_p}_{=p}
$$
so that 
$$
y_p=p-p''+p^{(4)}+...+(-1)^{m-1}p^{(2m-2)}
$$
For greater generality you can now add all the zero terms of the higher derivatives to get an infinite sum.

Or in short with the infinite geometric series
$$
\frac1{1+D^2}=\sum_{n=0}^\infty(-D^2)^n.
$$
In any application you have to care for the convergence, applied to a polynomial that is not a problem as the series is a finite sum.
