# Solving 4 coupled linear 1st order ODE with constant coefficients

I'm getting trouble solving these 4 coupled differential equations, just can't get a result while it keep running.

eqs =
{
A1'[t] ==  I (μ0 - 2*g) A1[t] - I g0 (B1[t] + C1[t]),
B1'[t] == -I (g0 (2 + n) A1[t] + Ω C1[t] + g0 D1[t]),
C1'[t] == -I (g0 (2 + n) A1[t] + Ω B1[t] + g0 D1[t]),
D1'[t] == -I g0 (1 + n) (B1[t] + C1[t]) - I (μ0 - 2 g) D1[t],
A1[0] == 0, B1[0] == 0, C1[0] == 0, D1[0] == 1
};
dusol = DSolve[eqs, {A1, B1, C1, D1}, t][[1]] // FullSimplify

• Don't use \!$$\*SuperscriptBox[\(A1$$, $$'$$]\)[t] to describe the derivative. Use A1'[t] instead. Sep 16, 2020 at 7:28
• Good question. Something definitely is wrong! Sep 17, 2020 at 0:15
• @xzczd The problem as stated runs forever for Version 12.1.1. Eliminating both the initial conditions and FullSimplify yields an enormous answer after several minutes. The matrix of the right sides of the equations is not Hermitian, but Eigensystem does produce results quickly. However, JordanDecomposition crashes. Possibly a bug. Sep 17, 2020 at 2:28
• …You got this warning with \!$$\*SuperscriptBox[\(A1$$, $$'$$]\)[t], right? I guess you obtained this by stroking Ctrl + 6 first? Then you've typed ' in wrong way. You should directly type ' i.e. A1'[t] is the correct input. Now copy the code in your question back to Mathematica and retry. (BTW, if you insist on stroking Ctrl + 6 first, then what you should type is \[Prime]. ) Sep 17, 2020 at 2:43
• Problem with JordanDecomposition noted in earlier comment report to Wolfram, Inc as CASE:4633110. Sep 17, 2020 at 13:50

I've found 3 work-arounds.

First one is to make use of the ConstantsGrouping  package:

head = {A1, B1, C1, D1};

{newrhs, rule} = GroupConstants[eqs[[;; 4, -1]], {t, var} // Flatten]

lhs = eqs[[;; 4, 1]]

solmid = DSolve[{lhs == newrhs // Thread, eqs[[5 ;;]]}, head, t]; // AbsoluteTiming
(* {12.4674, Null} *)

sol = var /. solmid[[1]] /. rule;


Second one is based on LaplaceTransform:

teqs = LaplaceTransform[eqs[[;; 4]], t, s] /. Rule @@@ eqs[[5 ;;]]

tsol = Solve[teqs, LaplaceTransform[var, t, s]][[1, All, -1]]

sol2 = InverseLaplaceTransform[tsol, s, t]; // AbsoluteTiming
(* {0.968162, Null} *)


Last, and the fastest one is to utilize MatrixExp:

{barray, marray} = CoefficientArrays[eqs[[;; 4, -1]], var];

sol3 = MatrixExp[marray t, eqs[[5 ;;, -1]]]; // AbsoluteTiming
(* {0.458083, Null} *)
`
• Thank you very much for your solutions. They're very useful. I've solved my problem. And importantly, I left out the initial condition $μ0=2g,g0=g$! Sep 19, 2020 at 7:51