Solution of $y''+xy=0$

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)$$

$$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$$

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.