# 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 -> c, y' -> 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}+O\left(x^9\right)$$

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 == 1, y == 1, u == 1, v == 1,
x' == 0, y' == 0, u' == 0, v' == 0},
{x, y, u, v}, {t, 0, 100}]


$$\left\{\left\{u\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),v\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),x\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),y\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right)\right\}\right\}$$

one can simplify to get trigonometric expression.

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


Then plotting yields,

 a = 0.2; k = 0.5;
GraphicsGrid[{
{Plot[{x[t] /. solution[]}, {t, 0, 10}],
Plot[{y[t] /. solution[]}, {t, 0,
10}]}, {Plot[{u[t] /. solution[]}, {t, 0, 10}],
Plot[{v[t] /. solution[]}, {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)

• Why not use AsymptoticDSolveValue? Jan 30 at 1:59
• @CarlWoll, AsymptoticDSolveValue is available from v12. Jan 30 at 12:08
• Actually, it is available starting in 11.3. Jan 30 at 23:56