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)