equation solving - General solution for a linear ODE set with complicated coefficient

This is the original problem that motivated me to ask this question. I encountered it when trying to reproduce the deduction in this paper. (I'll paste the relevant part below to make this question self-contained, of course.)

In page 4874, the author deduced a set of linear ODE

where $\Theta^*$, $U_x^*$, $U_y^*$, $U_z^*$ are unknown functions of $z$, $D^n()$ denotes $\frac{d^n}{dz^n}$ i.e. $D^n()\Theta^*$ evaluates to $\frac{d^n\Theta^*}{dz^n}$, $\kappa$, $p$, $q$, $w$, $h$, $\epsilon$, $m$ are constants.

Then the author found its general solution in a hard-to-understand way (at least for me), so I decided to try solving it with DSolve, but a direct solving seems to be way too slow:

coef = ( {

    {(k (p^2 q^2 + p^2 w^2 - p/h + d[2]))/m^2, e p^2 q/h, e p^2 w/h, e p/h d@1},
    {k p q, p^2 w^2 + p^2 m^2 (q^2 - 1) + d@2, p^2 q (m^2 - 1) w, p q (m^2 - 1) d@1},
    {k p w, p^2 q w (m^2 - 1), p^2 q^2 + p^2 m^2 (w^2 - 1) + d@2, p w (m^2 - 1) d@1},
    {k d@1, p q (m^2 - 1) d@1, p w (m^2 - 1) d@1, p^2 (q^2 + w^2 - m^2) + m^2 d@2}
   } );    
var = u[#][z] & /@ Range@4;
eqoriginal = Expand[coef.var == 0 // Thread] /. d[n_] u[i_][z] :> Derivative[n][u@i][z]

DSolve[eqoriginal, var, z]
(* Never finished *)

At this point it's natural to guess it's the complicated coefficients that slow down the solving, so I replaced them with simpler symbol, but this approach is still way too slow:

eq = Expand[(SparseArray[{i_, i_} -> d@2, {4, 4}, 0] Array[c1, {4, 4}] + 
        SparseArray[{{4, 4} -> 1, {_, 4} | {4, _} -> d@1}, {4, 4}, 1] Array[
          c2, {4, 4}]).var == 0 // Thread] /. u[i_][z] d[n_] -> Derivative[n][u[i]][z];

DSolve[eq, var, z] // AbsoluteTiming
(* Never finished *)

How can I solve the ODE set with Mathematica, in a relatively short period of time?

Answer

Because the equations are linear, homogeneous, and with constant coefficients, the solutions are sums of exponentials, possibly multiplied by powers of z, the independent variable. These terms can be obtained by constructing the determinant of the coefficients, converting it to a single ODE, and solving it.

Collect[Det[coef /. d[2] -> d[1]^2], d[1], Simplify[#] u[z] &] /. 
    d[1]^n_ u[z] -> D[u[z], {z, n}];
DSolve[% == 0, u[z], z][[1, 1]];
sols = Cases[%, a_ C[_] -> a, Infinity]
(* {E^((Sqrt[p/h + (e p)/h + p^2 - (p Sqrt[1 + 2 e + e^2 - 2 h p + 2 e h p + h^2 p^2])/h - 2 p^2 q^2 - 2 p^2 w^2] z)/Sqrt[2]), 
    E^(-((Sqrt[p/h + (e p)/h + p^2 - (p Sqrt[1 + 2 e + e^2 - 2 h p + 2 e h p + h^2 p^2])/h - 2 p^2 q^2 - 2 p^2 w^2] z)/Sqrt[2])), 
    E^((Sqrt[p/h + (e p)/h + p^2 + (p Sqrt[1 + 2 e + e^2 - 2 h p + 2 e h p + h^2 p^2])/h - 2 p^2 q^2 - 2 p^2 w^2] z)/Sqrt[2]), 
    E^(-((Sqrt[p/h + (e p)/h + p^2 + (p Sqrt[1 + 2 e + e^2 - 2 h p + 2 e h p + h^2 p^2])/h - 2 p^2 q^2 - 2 p^2 w^2] z)/Sqrt[2])), 

    E^(p Sqrt[m^2 - q^2 - w^2] z), E^(p Sqrt[m^2 - q^2 - w^2] z) z, 
    E^(-p Sqrt[m^2 - q^2 - w^2] z), E^(-p Sqrt[m^2 - q^2 - w^2] z) z} *)

Thus, we would expect the four dependent variables to be sums of these eight functions, multiplied by constants to be determined. Doing so, we find that the sixth and eighth elements of sols are spurious, in the sense that their coefficients are zero. ({c1, c2, c3, c4} are the coefficients for {u1, u2, u3, u4}, respectively.

Unevaluated[Expand[coef.{c1, c2, c3, c4} y[z]] /. {d[1] y[z] -> D[y[z], z], 
    d[2] y[z] -> D[y[z], {z, 2}]}] /. y[z] -> sols[[6]];
    op = Simplify[%]; Solve[Thread[op == 0], {c1, c2, c3, c4}] // Simplify
(* {{c1 -> 0, c2 -> 0, c3 -> 0, c4 -> 0}} *)

and similarly for sols[[8]]. Coefficients for the fifth element are

Unevaluated[Expand[coef.{c1, c2, c3, c4} y[z]] /. {d[1] y[z] -> D[y[z], z], 
    d[2] y[z] -> D[y[z], {z, 2}]}] /. y[z] -> sols[[5]];
    op = Simplify[%]; Solve[Thread[op == 0], {c1, c2, c3, c4}] // Simplify
(* {{c1 -> 0, c4 -> (-c2 q - c3 w)/Sqrt[m^2 - q^2 - w^2]}} *)

and for the seventh

(* {{c1 -> 0, c4 -> (c2 q + c3 w)/Sqrt[m^2 - q^2 - w^2]}} *)

In both of the cases, two of the four coefficients are arbitrary. Coefficients for the remaining four elements are more complicated, and we give only one here.

Unevaluated[Expand[coef.{c1, c2, c3, c4} y[z]] /. {d[1] y[z] -> D[y[z], z], 

    d[2] y[z] -> D[y[z], {z, 2}]}] /. y[z] -> sols[[1]];
op = Simplify[%]; Solve[Thread[op[[2 ;; 4]] == 0], {c2, c3, c4}] // Simplify
(* {{c2 -> (2 c1 h k q)/(m^2 (-1 - e + h p + Sqrt[e^2 + (-1 + h p)^2 + 2 e (1 + h p)])), 
     c3 -> (2 c1 h k w)/(m^2 (-1 - e + h p + Sqrt[e^2 + (-1 + h p)^2 + 2 e (1 + h p)])), 
     c4 -> (Sqrt[2] c1 h k Sqrt[-((p (-1 - e - h p + Sqrt[e^2 + (-1 + h p)^2 + 2 e (1 + h p)] + 2 h p q^2 + 2 h p w^2))/h)])/(m^2 p (-1 - e + h p + Sqrt[e^2 + (-1 + h p)^2 + 2 e (1 + h p)]))}} *)

Thus, one of the four coefficients is arbitrary, as expected. In all, therefore, the solution consists of six functions (elements 1 - 5, 7 of sols) and eight arbitrary constants.

Obtaining this solution required negligible computer time. The same cannot be said of my time.

This system of ODEs also can be solved by converting them to a single ODE, which is only sixth-order, due to cancellations. Moreover, the operator factors into four-order and second-order operators, the former yielding elements 1 - 4 of sols, and the second elements 5 and 7.

Solution by Laplace Transform

Alternatively, the system can be solved by Laplace transform, in which the determinant of the transformed system, equivalent to Det[coef], is a cubic in s^2, and the initial conditions are linear combinations the constants of integration. In brief, the derivation is as follows.

var = {u1[z], u2[z], u3[z], u4[z]};
lvar = Thread[LaplaceTransform[var, z, s]];
Thread[Simplify[Expand[coef.var] /. 
    {d[1] u_[z] -> D[u[z], z], d[2] u_[z] -> D[u[z], {z, 2}]}] == 0];
LaplaceTransform[%, z, s];
{ls, lcoef} = CoefficientArrays[%, lvar] // Normal // Simplify
(* {{-((e m^2 p u4[0] + h k (s u1[0] + Derivative[1][u1][0]))/(h m^2)), 
     -s u2[0] - (-1 + m^2) p q u4[0] - Derivative[1][u2][0], 
     -s u3[0] - (-1 + m^2) p w u4[0] - Derivative[1][u3][0], 

     -k u1[0] - (-1 + m^2) p (q u2[0] + w u3[0]) - m^2 (s u4[0] + Derivative[1][u4][0])}, 
    {{(k (-p + h s^2 + h p^2 (q^2 + w^2)))/(h m^2), (e p^2 q)/h, (e p^2 w)/h, (e p s)/h}, 
     {k p q, s^2 + p^2 (m^2 (-1 + q^2) + w^2), (-1 + m^2) p^2 q w, (-1 + m^2) p q s}, 
     {k p w, (-1 + m^2) p^2 q w, p^2 q^2 + s^2 + m^2 p^2 (-1 + w^2), (-1 + m^2) p s w}, 
     {k s, (-1 + m^2) p q s, (-1 + m^2) p s w, m^2 s^2 + p^2 (-m^2 + q^2 + w^2)}}} *)
inv = Inverse[lcoef] // Simplify;

As expected, lcoef is coef with d[1] and d[2] replaced by s and s^2, respectively, and ls is a vector of initial conditions. The LaplaceTransform of the solution is -inv.ls. Poles in this quantity are the arguments of the exponentials in sols, defined previously. Consistent with the earlier derivation, the LaplaceTransform of u1 (only) does not have poles for elements 5 and 7 of sols. The actual solution,

InverseLaplaceTransform[-inv.ls, s, z] // Simplify

with a LeafCount of 13057, is too large to be reproduced here. Computation time totals about 150 sec.