# Complicated System of Second Order Differential Equations

I am trying to solve a coupled system of second order differential equations using either DSolve[] or NDSolve[]. I believe Mathematica may not be powerful enough to find the solution because when I use the same syntax and simpler equations, it solves them correctly but when I run the full code below, it gives me my input back.

Here is my code:

SetAttributes[epsilon, {Constant, Protected}]
SetAttributes[alpha, {Constant, Protected}]
SetAttributes[beta, {Constant, Protected}]

eqns = {phi''[t] ==
-(3 (a[t]^2 a'[t] phi'[t] +
epsilon a[t]^2 a'[t] phi'[t] +
epsilon a[t]^2 Abs[alpha]^2 a'[t] phi'[t] +
epsilon Abs[alpha]^2 b[t]^2 b'[t] phi'[t]))/(a[t]^3 + epsilon a[t]^3 +
epsilon a[t]^3 Abs[alpha]^2 + epsilon Abs[alpha]^2 b[t]^3),
chi''[t] ==
-(3 (epsilon a[t]^2 Abs[beta]^2 a'[t] chi'[t] +
b[t]^2 b'[t] chi'[t] + epsilon b[t]^2 b'[t] chi'[t] +
epsilon Abs[beta]^2 b[t]^2 b'[t] chi'[t]))/(epsilon a[t]^3 Abs[beta]^2 +
b[t]^3 + epsilon b[t]^3 + epsilon Abs[beta]^2 b[t]^3),
a''[t] ==
1/(4 a[t]) (-2a'[t]^2 + epsilon a[t]^2 Abs[beta]^2 chi'[t]^2 +
a[t]^2 phi'[t]^2 + epsilon a[t]^2 phi'[t]^2 +
epsilon a[t]^2 Abs[alpha]^2 phi'[t]^2),
b''[t] ==
1/(4 b[t]) (-2 b'[t]^2 + b[t]^2 chi'[t]^2 +
epsilon b[t]^2 chi'[t]^2 + epsilon Abs[beta]^2 b[t]^2 chi'[t]^2 +
epsilon Abs[alpha]^2 b[t]^2 phi'[t]^2)};

sol = DSolve[eqns, {phi[t], chi[t], a[t], b[t]}, t]


I've also tried putting in initial conditions and get the same result. Is it possible Mathematica just can't handle a system this complicated? Is there anything I can do to make it work? Thanks!

• You have a syntax error: "beta]^2[a]'[t] " What should 2[a] mean? Aug 18, 2022 at 19:19
• Oops, thank you. It's fixed now. Aug 18, 2022 at 19:22
• For NDSolve you need initial conditions and numeric values for the parameters. For DSolve it would not be surprising to me if a nonlinear, high-dimensional system cannot be solved symbolically (by DSolve). Aug 18, 2022 at 19:27

Here's a start with NDSolve. As Michael stated, need to assign values to constants and set up some initial conditions. You can change them as you like although I had to adjust them a bit before I was able to successfully integrate from t=0 to t=2 without encountering a stiff system suggesting a singularity was reached such as division by zero:

epsilon = 1;
alpha = 1.5;
beta = -2.25;

initVals = {phi[0] == 1/10, phi'[0] == 1/4, chi[0] == 1, chi'[0] == 1,
a[0] == 1/4, a'[0] == 1/4, b[0] == -1, b'[0] == -1};
eqns = {phi''[
t] == -(3 (a[t]^2 a'[t] phi'[t] + epsilon a[t]^2 a'[t] phi'[t] +
epsilon a[t]^2 Abs[alpha]^2 a'[t] phi'[t] +
epsilon Abs[alpha]^2 b[t]^2 b'[t] phi'[t]))/(a[t]^3 +
epsilon a[t]^3 + epsilon a[t]^3 Abs[alpha]^2 +
epsilon Abs[alpha]^2 b[t]^3),
chi''[t] == -(3 (epsilon a[t]^2 Abs[beta]^2 a'[t] chi'[t] +
b[t]^2 b'[t] chi'[t] + epsilon b[t]^2 b'[t] chi'[t] +
epsilon Abs[beta]^2 b[t]^2 b'[t] chi'[t]))/(epsilon a[
t]^3 Abs[beta]^2 + b[t]^3 + epsilon b[t]^3 +
epsilon Abs[beta]^2 b[t]^3),
a''[t] ==
1/(4 a[t]) (-2 a'[t]^2 + epsilon a[t]^2 Abs[beta]^2 chi'[t]^2 +
a[t]^2 phi'[t]^2 + epsilon a[t]^2 phi'[t]^2 +
epsilon a[t]^2 Abs[alpha]^2 phi'[t]^2),
b''[t] ==
1/(4 b[t]) (-2 b'[t]^2 + b[t]^2 chi'[t]^2 +
epsilon b[t]^2 chi'[t]^2 +
epsilon Abs[beta]^2 b[t]^2 chi'[t]^2 +
epsilon Abs[alpha]^2 b[t]^2 phi'[t]^2)};
{f1, f2, f3, f4} =
NDSolveValue[Join[{eqns, initVals}], {phi, chi, a, b}, {t, 0, 2}]
Plot[{f1[t], f2[t], f3[t], f4[t]}, {t, 0, 2}]