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)