How to make a code to find Taylor series symbolic solution to four coupled nonlinear differential equations?

I am trying to modify the existing code developed by Michael E2 in this question here. His solution was for one differential equation. I like his code because it has ability to solve nonlinear differential equation symbolically. I would like to extend it to be able to solve four coupled differential equation of second order.
I am running the following test codes and my attempt to modify the existing code:

ClearAll[x, t, a, b, c, xx];
(*******       testing for simple pendulum equation   ***************)
(*******           Series for one SONDE                 *************)
Clear[seriesDSolve];
seriesDSolve[ode_,y_,iter:{x_,x0_,n_},ics_:    {}]:=Module[{dorder,coeff},dorder=Max[0,Cases[ode,Derivative[m_][y][x]:>m,Infinity]];
(coeff=Fold[Flatten[{#1,Solve[#2==0/.#1,Union@Cases[#2/.#1,Derivative[m_][y][x0]/;m>=dorder,Infinity]]}]&,ics,Table[SeriesCoefficient[ode/.Equal->Subtract,{x,x0,k}],{k,0,Max[n-   dorder,0]}]];
Series[y[x],iter]/.coeff)/;dorder>0]

seriesDSolve[
y''[x] + a*Sin[y[x] == 0], y, {x, 0, 8}, {y[0] -> c, y'[0] -> 0}

and the output is

$c-frac{1}{2} x^2 (a sin (c))+frac{1}{24} a^2 x^4 sin (c) cos (c)+frac{1}{720} x^6 left(3 a^3 sin ^3(c)-a^3 sin (c) cos ^2(c)right)+frac{a x^8 left(a^3 sin (c) cos ^3(c)-33 a^3 sin ^3(c) cos (c)right)}{40320}+Oleft(x^9right)$

which is correct.

Now a test code for 4 coupled differential equations (it should also handle nonlinear ones as well).

 ClearAll[k,a];
 solution = DSolve[{x''[t] + a*x[t] + k*v[t] == 0,
                                y''[t] + a*v[t] + k*x[t] == 0, 
  u''[t] + a*v[t] + k*x[t] == 0,
  v''[t] + a*v[t] + k*x[t] == 0,
  x[0] == 1, y[0] == 1, u[0] == 1, v[0] == 1,
  x'[0] == 0, y'[0] == 0, u'[0] == 0, v'[0] == 0},
  {x, y, u, v}, {t, 0, 100}]

and the answer is :

$left{left{uto left({t} {f4a1}frac{1}{2} e^{t left(-sqrt{-a-k}right)} left(e^{2 t sqrt{-a-k}}+1right)right),vto left({t} {f4a1}frac{1}{2} e^{t left(-sqrt{-a-k}right)} left(e^{2 t sqrt{-a-k}}+1right)right),xto left({t} {f4a1}frac{1}{2} e^{t left(-sqrt{-a-k}right)} left(e^{2 t sqrt{-a-k}}+1right)right),yto left({t} {f4a1}frac{1}{2} e^{t left(-sqrt{-a-k}right)} left(e^{2 t sqrt{-a-k}}+1right)right)right}right}$

one can simplify to get trigonometric expression.

 p = FullSimplify[{x[t], y[t], u[t], v[t]} /. solution[[1]]] 

Then plotting yields,

 a = 0.2; k = 0.5;
 GraphicsGrid[{
{Plot[{x[t] /. solution[[1]]}, {t, 0, 10}], 
 Plot[{y[t] /. solution[[1]]}, {t, 0, 
 10}]}, {Plot[{u[t] /. solution[[1]]}, {t, 0, 10}], 
 Plot[{v[t] /. solution[[1]]}, {t, 0, 10}]}}]

Then working on series for 4 coupled SODE’s,

Clear[seriesDSolve];
seriesDSolve[ode_, f_, iter : {t_, t0_, n_}, ics_: {}] := 
Module[{dorder, coeff}, 
dorder = Max[0, Cases[ode, Derivative[m_][f][t] :> m, Infinity]];
(coeff = 
 Fold[Flatten[{#1, 
     Solve[#2 == 0 /. #1, 
      Union@Cases[#2 /. #1, Derivative[m_][f][t0] /; m >= dorder, 
        Infinity]]}] &, ics, 
  Table[SeriesCoefficient[
    ode /. Equal -> Subtract, {t, t0, k}], {k, 0, 
    Max[n - dorder, 0]}]];
  Series[f[t], iter] /. coeff) /; dorder > 0]

Unfortunately, the modified code does not work. My list approach is not correct.
Any suggestions on what I am missing?
(MMA v11)