Solution of $y''+xy=0$ The differential equation $y''+xy=0$ is given.
Find the solution of the differential equation, using the power series method.
That's what I have tried:
We are looking for a solution of the form $y(x)= \sum_{n=0}^{\infty} a_n x^n$ with radius of convergence of the power series $R>0$.
Then:
$$y'(x)= \sum_{n=1}^{\infty} n a_n x^{n-1}= \sum_{n=0}^{\infty} (n+1) a_{n+1} x^n$$
$$y''(x)= \sum_{n=1}^{\infty} (n+1) n a_{n+1} x^{n-1}= \sum_{n=0}^{\infty} (n+2) (n+1) a_{n+2} x^n$$
Thus:
$$\sum_{n=0}^{\infty} (n+2) (n+1) a_{n+2} x^n+ x \sum_{n=0}^{\infty} a_n x^n=0 \\ \Rightarrow \sum_{n=0}^{\infty} (n+2)(n+1) a_{n+2} x^n+ \sum_{n=0}^{\infty} a_n x^{n+1}=0 \\ \Rightarrow \sum_{n=0}^{\infty} (n+2)(n+1) a_{n+2} x^n+ \sum_{n=1}^{\infty} a_{n-1} x^n=0 \\ \Rightarrow 2a_2+\sum_{n=1}^{\infty} \left[ (n+2) (n+1) a_{n+2}+ a_{n-1}\right] x^n=0$$
So it has to hold:
$$a_2=0 \\ (n+2) (n+1) a_{n+2}+a_{n-1}=0, \forall n=1,2,3, \dots$$
For $n=1$: $3 \cdot 2 \cdot a_3+ a_0=0 \Rightarrow a_3=-\frac{a_0}{6}$
For $n=2$: $4 \cdot 3 \cdot a_4+a_1=0 \Rightarrow a_4=-\frac{a_1}{12}$
For $n=3$: $5 \cdot 4 \cdot a_5+a_2=0 \Rightarrow a_5=0$
For $n=4$: $6 \cdot 5 \cdot a_6+a_3=0 \Rightarrow 30 a_6-\frac{a_0}{6}=0 \Rightarrow a_6=\frac{a_0}{6 \cdot 30}=\frac{a_0}{180}$
For $n=5$: $7 \cdot 6 \cdot a_7+ a_4=0 \Rightarrow 7 \cdot 6 \cdot a_7-\frac{a_1}{12}=0 \Rightarrow a_7=\frac{a_1}{12 \cdot 42}$
Is it right so far? If so, how could we find a general formula for the coefficients $a_n$?
EDIT:  Will it be as follows:
$$a_{3k+2}=0$$
$$a_{3k}=(-1)^k \frac{a_0}{(3k)!} \prod_{i=0}^{k-1} (3i+1)$$
$$a_{3k+1}=(-1)^k \frac{a_1}{(3k+1)!} \prod_{i=0}^{k-1} (3i+2)$$
If so, then do we have to write seperately the formula for the coefficients of $x^0, x^1$, because otherwhise the sum would be from $0$ to $-1$ ?
 A: Follow-up from my comment, every $a_{n+3}$ can be expressed in terms of $a_n$.
$$a_{n+3} = -\frac{a_{n}}{(n+3)(n+2)}$$ 
Since $a_2 = 0$, all coefficients with $n = 3k + 2$ will be $0$
The rest will have to based on the initial conditions of $a_0 = y(0)$ and $a_1 = y'(0)$
For $n = 3k$
$$a_3 = -\frac{a_0}{3\cdot2}$$
$$a_6 = -\frac{a_3}{6\cdot5} = \frac{a_0}{6\cdot5\cdot3\cdot2} $$
This can be extrapolated to
$$ a_{n} = \pm\frac{a_0}{n(n-1)(n-3)(n-4)\cdots\cdot6\cdot5\cdot3\cdot2} $$
Basically, the denominator is a product of all integers from $1$ to $n$, skipping every 3rd number and alternating signs. An alternate notation is
$$a_{3k} = (-1)^k\frac{a_0}{(3k)!}\prod_{i=0}^{k-1} (3i+1) $$
The same can be done for $n = 3k+1$
$$a_{3k+1} = (-1)^k\frac{a_1}{(3k+1)!}\prod_{i=0}^{k-1} (3i+2)$$
A: $$y''+xy=0$$
$$y=\sum_{n=0}^{\infty} c_n x^n, y'=\sum_{n=1}^{\infty} nc_n x^{n-1}, y''=\sum_{n=2}^{\infty} n(n-1)c_n x^{n-2}$$
$$\therefore y''+xy=0=\underbrace{\sum_{n=2}^{\infty} n(n-1)c_n x^{n-2}}_{k=n-2\Rightarrow n=k+2}+\underbrace{\sum_{n=0}^{\infty} c_n x^{n+1}}_{k=n+1\Rightarrow n=k-1}=\sum_{k=0}^{\infty} (k+2)(k+1)c_{k+2} x^{k}+\sum_{k=1}^{\infty} c_{k-1} x^k$$
$$=2c_2+\sum_{k=1}^{\infty} (k+2)(k+1)c_{k+2} x^{k}+\sum_{k=1}^{\infty} c_{k-1} x^k=2c_2+\sum_{k=1}^{\infty} [(k+2)(k+1)c_{k+2}+ c_{k-1}] x^k=0$$
Thus,
$$2c_2=0\Rightarrow c_2=0$$
$$(k+2)(k+1)c_{k+2}+ c_{k-1}=0\Rightarrow c_{k+2}=-\frac{c_{k-1}}{(k+2)(k+1)}\forall k=1,2,3,\ldots$$
Choosing $c_0=1$ and $c_1=0$, we find
$$c_2=0,c_3=-1/6,c_4=0, c_5=0, c_6=1/180\ldots$$
and so on.
Choosing $c_0=0$ and $c_1=1$, we find
$$c_2=0,c_3=0,c_4=-1/12, c_5=0, c_6=0,c_7=1/504\ldots$$
and so on.
Thus, the two solutions are
$$y_1=1-\frac{1}{6}x^3+\frac{1}{180}x^6+\ldots$$
$$y_2=x-\frac{1}{12}x^4+\frac{1}{504}x^7+\dots$$
A: Disclaimer. This is not an answer, but rather a too long comment with some graphics in it.
The equation looks like the common ODE for a harmonic oscillator: $y'' + \omega^2 y = 0$ with the square of the frequency varying proportional to "time" $x$ . Numerical simulation with the initial conditions $y(0)=0$ and $y'(0)=1$ reveals that the solution indeed looks like that (I hate solutions without a picture, you know): 
The viewport is $\,0 < x < 40\,$ and $\,-1.5 < y < +1.5\,$.
Program (Delphi Pascal) snippet for doing the calculations and the drawing (not optimized at all):

{
The equations of motion are solved numerically as follows.
Start with: y(0) = 0 ; x = 0 ;
  (y(dx) - y(0))/dx = v  ==>  y(dx) = y(0) + v.dx
  y(x + dx) - 2.y(x) + y(x - dx)
  ------------------------------ + x.y(x) = 0   ==>
               dx^2
  y(x + dx) = 2.y(x) - y(x - dx) - dx^2.x.y(x)  ==>
  y(0.dx) = 0
  y(1.dx) = y(0.dx) + v.dx
  y(2.dx) = 2.y(1.dx) - y(0.dx) - dx^2.x.y(1.dx)
  y(3.dx) = 2.y(2.dx) - y(1.dx) - dx^2.x.y(2.dx)
  ..............................................
  y(k+1).dx) = 2.y(k.dx) - y((k-1).dx) - dx^2.x.y(k.dx)
}
procedure bereken;
const
  N : integer = 40000;
  v : double = 1;
var
  x,dx,y0,y1,y2 : double;
  k : integer;
begin
  x := 0; dx := 0.001;
  y0 := 0; y1 := y0+v*dx;
  Form1.Image1.Canvas.MoveTo(x2i(0),y2j(0));
  for k := 0 to N-1 do
  begin
    x := x + dx;
    y2 := 2*y1 - y0 - dx*dx*x*y1;
    Form1.Image1.Canvas.LineTo(x2i(x),y2j(y2));
    y0 := y1; y1 := y2;
  end;
end;

Note that the solution becomes very oscillatory (i.e. singular) for $x\to\infty$ . However, the frequency only varies with the square root of the distance: $\omega=\sqrt{x}$ , therefore it doesn't happen immediately.
