To compute a power series up to any desired order, you can use the tools of Taylor series arithmetic that are a part of the Automatic/Algorithmic Differentiation theory. As differentiation of a power series contains an index shift, new coefficients result on the left side while evaluating or updating the evaluation on the right side.
To simplify the Taylor series arithmetic introduce $u=\cos y$, $v=\sin y$. Then
\begin{align}
u'&=-vy'\\
v'&=~~~uy'\\
y''&=-axu.
\end{align}
This is now a problem where the multiplication of Taylor series is the most complicated operation which allows to establish the equations for the coefficients of these 3 series
$$y=\sum_{n=0}^\infty y_kx^k, ~~~u=\sum_{n=0}^\infty u_kx^k,~~~ v=\sum_{n=0}^\infty v_kx^k.
$$
These can then be computed by successively inserting the already computed coefficients,
\begin{align}
y_0&=y(0)& u_0&=\cos(y_0),& v_0&=\sin(y_0),\\
y_1&=y'(0)& u_1&=-v_0y_1,& v_1&=u_0y_1,\\
2y_2&=0,& 2u_2&=-v_0(2y_2)-v_1y_1,& 2v_2&=u_0(2y_2)+u_1y_1,\\
6y_3&=-au_0,& 3u_3&=-\sum_{k=1}^3v_{3-k}(ky_k),& 3v_3&=\sum_{k=1}^3u_{3-k}(ky_k),\\
&\vdots\\
n(n-1)y_n&=-au_{n-3},&nu_n&=-\sum_{k=1}^nv_{n-k}(ky_k),& nv_n&=\sum_{k=1}^nu_{n-k}(ky_k),\\
\end{align}
For further information on the automatic creation of such algorithms, see
- FADBAD++ by Stauning and Bendtsen (project active 1996-2003), the third example in the crash-course section of their description site,
- VNODE, used for examples in papers by N. Nedialkov (2001)
- TSM - ODE playground by user m4r35n357 (active 2018-), Taylor series manipulation with a focus on ODE