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.
