Initial value differential equation solution using series method about $x=0$ 
If $y''-15xy'+8y=e^{x}$ and $y(0)=y'(0)=-2$. Then find first $5$ non zero terms of  series solution of that equation about $x=0$

What i try :: let $\displaystyle y=\sum^{\infty}_{n=0}b_{n}x^n$ and $\displaystyle y'=\sum^{\infty}_{n=1}nb_{n}x^{n-1}$ and $\displaystyle y''=\sum^{\infty}_{n=2}n(n-1)b_{n}x^{n-2}$ and $b_{0}=b_{1}=-2$.
Then equation convert into
$$\sum^{\infty}_{n=2}n(n-1)b_{n}x^{n-2}-15\sum^{\infty}_{n=1}nb_{n}x^{n}+8\sum^{\infty}_{n=0}b_{n}x^{n}=e^{x}$$
$$\sum^{\infty}_{n=0}(n+1)(n+2)b_{n+2}x^{n}-15\sum^{\infty}_{n=1}nb_{n}x^{n}+8\sum^{\infty}_{n=0}b_{n}x^{n}=e^{x}$$
$$b_{0}+2b_{2}+\sum^{\infty}_{n=1}\bigg[(n+1)(n+2)b_{n+2}-15nb_{n}+8b_{n}\bigg]=e^{x}$$
Now i did not understand How can i proceed further.
help me please. Thanks
 A: Suppose that you have a series
$$y(x) = \sum_{n=0}^{+\infty} b_nx^n$$
with a positive radius of convergence, which is solution of the equation. Then in the open disc of convergence, you can differentiate term by term, so you get
$$y''(x) - 15xy'(x) + 8y(x) = \sum_{n=2}^{+\infty} n(n-1)b_nx^{n-2} - 15x \sum_{n=1}^{+\infty}nb_nx^{n-1} + 8 \sum_{n=0}^{+\infty} b_n x^n$$
$$=\sum_{n=0}^{+\infty} (n+1)(n+2)b_{n+2}x^{n} - 15 \sum_{n=1}^{+\infty}nb_nx^{n} + 8 \sum_{n=0}^{+\infty} b_n x^n $$
$$= 2b_2 + 8b_0 + \sum_{n=1}^{+\infty} \left[ (n+1)(n+2)b_{n+2}+(8-15n)b_n\right]x^n$$
This has to be equal to $e^x$, i.e. you have
$$2b_2 + 8b_0 + \sum_{n=1}^{+\infty} \left[ (n+1)(n+2)b_{n+2}+(8-15n)b_n\right]x^n = 1 + \sum_{n=1}^{+\infty} \frac{x^n}{n!}$$
Identifying the terms of the same degree, you get that
$$2b_2 + 8b_0 = 1 \quad \text{and, for all }n\geq 1, (n+1)(n+2)b_{n+2}+(8-15n)b_n = \frac{1}{n!}$$
Now because $y(0)=y'(0)=-2$, then $b_0=b_1=-2$, so you deduce
$$b_2 = \frac{17}{2}$$
then
$$b_3 = \frac{-13}{6}$$
then
$$b_4  = \frac{125}{8}$$
and finally
$$b_5  = -4$$
A: A different method, that can be easily applied even if you don't know the series development of the RHS (but calculating its derivatives is more or less equivalent).
Suppose $y$ can be developed Maclaurin series
$$
y(x)=\sum_{k=0}^\infty\frac{y^{(k)}(0)}{k!}x^k
$$
From the equation itself we get
$$
y''=15xy'-8y+e^x\tag1
$$
so that
$$
y''(0)=15\cdot0\cdot(-2)-8\cdot(-2)+e^0=17
$$
Differentiating $(1)$
\begin{align}
y''' &= 15xy''+15y'-8y'+e^x=\\
  &= 15xy''+7y'+e^x\tag2
\end{align}
so that
$$
y'''(0)=15\cdot0\cdot(17)+7\cdot(-2)+e^0=-13
$$
Next, differentiate $(2),$ obtaining
\begin{align}
y^{(4)} &= 15xy'''+15y''+7y''+e^x=\\
  &=15xy'''+22y''+e^x
\end{align}
and
$$
y^{(4)}(0)=0+22\cdot17+1=375
$$
and so on...
Don't forget to divide each coefficient by $k!$
Eventually, the first $5$ terms are
$$
-2,\qquad-2,\qquad\frac{17}{2},\qquad-\frac{13}{6},\qquad\frac{375}{4!}=\frac{125}{8}
$$
A: You started well and reached to
$$
\sum\limits_{0\, \le \,n} {n\left( {n - 1} \right)b_{\,n} \,x^{\,n - 2} }  - 15\sum\limits_{0\, \le \,n} {n\,b_{\,n} \,x^{\,n} }  + 8\sum\limits_{0\, \le \,n} {\,b_{\,n} \,x^{\,n} }  = e^{\,x} 
$$
note that you have better and keep the sum to start from $0$: some of the initial terms are intrinsically null.
What you have to do now is to "realign" the sums in such a way that all of them shows $x^n$, in this way
$$
\eqalign{
  & \sum\limits_{0\, \le \,n - 2} {n\left( {n - 1} \right)b_{\,n} \,x^{\,n - 2} }  - 15\sum\limits_{0\, \le \,n} {n\,b_{\,n} \,x^{\,n} }  + 8\sum\limits_{0\, \le \,n} {\,b_{\,n} \,x^{\,n} }   \cr 
  & \sum\limits_{0\, \le \,n } {\left( {n + 2} \right)\left( {n + 1} \right)b_{\,n + 2} \,x^{\,n} }  - 15\sum\limits_{0\, \le \,n} {n\,b_{\,n} \,x^{\,n} }  + 8\sum\limits_{0\, \le \,n} {\,b_{\,n} \,x^{\,n} }  \cr} 
$$
where in the first sum we changed the index, while in the second we added the term $0\cdot \, b_0 \cdot x^0$ which is null.
And now you can start and equate each term with that of the expansion of $e^x$, taking in consideration that you know that $b_0 = b_1 = -2$
That is
$$
\left( {n + 2} \right)\left( {n + 1} \right)b_{\,n + 2}  + \left( {8 - 15n} \right)b_{\,n}  = {1 \over {n!}}\quad \left| {\,b_{\,1}  = b_{\,0}  =  - 2} \right.
$$
or
$$
b_{\,n + 2}  = {{\left( {15n - 8} \right)} \over {\left( {n + 2} \right)\left( {n + 1} \right)}}b_{\,n}  + {1 \over {\left( {n + 2} \right)!}}\quad \left| {\,b_{\,1}  = b_{\,0}  =  - 2} \right.
$$
or
$$
\left\{ \matrix{
  b_{\,0}  =  - 2 \hfill \cr 
  b_{\,1}  =  - 2 \hfill \cr 
  b_2  =  - {8 \over 2}\left( { - 2} \right) + {1 \over 2} = {{17} \over 2} \hfill \cr 
  b_{\,3}  = {7 \over 6}\left( { - 2} \right) + {1 \over 6} =  - {{13} \over 6} \hfill \cr 
  \quad \quad  \vdots  \hfill \cr}  \right.
$$
A: This is a pure mechanical procedure so I will let a script in order to obtain the desired results.
First we know that $e^x = \sum_{k=0}^{\infty}\frac{x^k}{k!}$ so defining $Y = \sum_{k=0}^n a_k x^k$ we develop the approximate relationship
$$
\cases{
Y''-15x Y' + 8Y = \sum_{k=0}^n\frac{x^k}{k!}\ \ \ (*)\\
Y(0) = -2\\
Y'(0) = -2}
$$
grouping for powers of $x$ for $n=5$ we get the relations
$$
\left\{
\begin{array}{rcl}
 8 a_0+2 a_2-1&=&0 \\
 -7 a_1+6 a_3-1&=&0 \\
 -22 a_2+12 a_4-\frac{1}{2}&=&0 \\
 -37 a_3+20 a_5-\frac{1}{6}&=&0 \\
 a_0&=&-2 \\
 a_1&=&-2 \\
\end{array}
\right.
$$
Note that from the $n$ coefficient relations obtainable from $(*)$ we follow only with the lower $n-2$ due to the presence of the two initial conditions. So at this point, after solving the linear system we obtain
$$
Y_5 = -4 x^5+\frac{125 x^4}{8}-\frac{13 x^3}{6}+\frac{17 x^2}{2}-2 x-2
$$
Follows the MATHEMATICA script.
n = 5;
ex = Sum[x^k/k!, {k, 0, n}];
Y = Sum[Subscript[a, k] x^k, {k, 0, n}];
A = Table[Subscript[a, k], {k, 0, n}];
res = D[Y, {x, 2}] - 15 x D[Y, x] + 8 Y - ex;
coefs = Take[CoefficientList[res, x], {1, n - 1}];
cond1 = (Y /. {x -> 0}) == -2;
cond2 = (D[Y, x] /. {x -> 0}) == -2;
conds = {cond1, cond2};
equs = Thread[coefs == 0];
equstot = Join[conds, equs]
sol = Solve[equstot, A][[1]]
Yx = Y /. sol

and to evaluate the accuracy
xmax = 1;
soly = NDSolve[{y''[x] - 15 x y'[x] + 8 y[x] == Exp[x], y[0] == -2, y'[0] == -2}, y, {x, 0, xmax}][[1]];
gr1 = Plot[Yx, {x, 0, xmax}, PlotStyle -> Dashed];
gr2 = Plot[Evaluate[y[x] /. soly], {x, 0, xmax}];
Show[gr1, gr2]

