One thing you can do is replace the exponential terms E^(-x h[i]) with a variable e[i], solve the system, and replace the exponential terms at the end. Taking away the exponentials limits what Mathematica can do with the system. If it can solve this more generic system, we still win. And since there is less that can be tried, it is likely to be faster.
Update: Polynomial equations are also easier to solve. Originally, I thought rational functions might be converted to polynomial equations, but apparently not. We can multiply eq3 through by the denominators, expand, and simplify. Since I’m changing eq3, here is the complete code.
Clear[a, b, c, d, h, y, z, l];
l = 4;
vars = Flatten[Table[{a[i], b[i], c[i], d[i]}, {i, l}]];
eq1 = {{x (x a[1] + x b[1] + d[1] (1 - 2 z[1]) +
c[1] (-1 + 2 z[1]))}, {x (x a[1] - x b[1] + 2 c[1] z[1] +
2 d[1] z[1])}} == {{1}, {0}};
eq2 = {{a[4]}, {c[4]}} == {{0}, {0}};
eq3 = Table[
y[1 + i] y[i] e[i] * (* multiply by common denominator *)
{{E^(x h[i]) x a[i] + E^(-x h[i]) x b[i] +
E^(-x h[i]) d[i] (1 + x h[i] - 2 z[i]) +
E^(x h[i]) c[i] (-1 + x h[i] + 2 z[i])}, {E^(x h[i]) x a[i] -
E^(-x h[i]) x b[i] + E^(-x h[i]) d[i] (-x h[i] + 2 z[i]) +
E^(x h[i]) c[
i] (x h[i] + 2 z[i])}, {-((E^(x h[i]) x a[i] (1 + z[i]))/
y[i]) + (E^(-x h[i]) x b[i] (1 + z[i]))/
y[i] + (E^(x h[i]) c[i] (2 - x h[i] - 4 z[i]) (1 + z[i]))/
y[i] + (E^(-x h[i]) d[i] (2 + x h[i] - 4 z[i]) (1 + z[i]))/
y[i]}, {-((E^(x h[i]) x a[i] (1 + z[i]))/
y[i]) - (E^(-x h[i]) x b[i] (1 + z[i]))/
y[i] + (E^(-x h[i]) d[i] (1 - x h[i]) (1 + z[i]))/
y[i] - (E^(x h[i]) c[i] (1 + x h[i]) (1 + z[i]))/y[i]}} ==
y[1 + i] y[i] e[i] * (* multiply by common denominator *)
{{x a[1 + i] + x b[1 + i] + d[1 + i] (1 - 2 z[1 + i]) +
c[1 + i] (-1 + 2 z[1 + i])}, {x a[1 + i] - x b[1 + i] +
2 c[1 + i] z[1 + i] +
2 d[1 + i] z[
1 + i]}, {-((x a[1 + i] (1 + z[1 + i]))/
y[1 + i]) + (x b[1 + i] (1 + z[1 + i]))/
y[1 + i] + (c[1 + i] (2 - 4 z[1 + i]) (1 + z[1 + i]))/
y[1 + i] + (d[1 + i] (2 - 4 z[1 + i]) (1 + z[1 + i]))/
y[1 + i]}, {-((x a[1 + i] (1 + z[1 + i]))/
y[1 + i]) - (x b[1 + i] (1 + z[1 + i]))/
y[1 + i] - (c[1 + i] (1 + z[1 + i]))/
y[1 + i] + (d[1 + i] (1 + z[1 + i]))/y[1 + i]}}, {i, l - 1}];
(soln = Solve[
Join[{eq1, eq2},
eq3 /. {E^(-x h[i_]) :> 1/e[i], E^(x h[i_]) :> e[i]} //
Expand // Simplify],
vars]) // Short // AbsoluteTiming
The solution is quite big and chokes Simplify:
(solnExp = soln /. e[i_] :> E^(x h[i])) // ByteCount // AbsoluteTiming
(*
{15.650537, 729043792}
*)
It’s interesting that it takes longer to back substitute, than to solve.
We can simplify piecemeal, working from the bottom up. We can map Simplify to the parts at some level from the bottom (a sort of divide-and-conquer strategy). I use Fold to do it sequentially for levels -6, -8, -10, -12, printing out the timings and bytecounts for each level. The size is quickly reduced.
Depth[solnExp]
(*
38
*)
soln1 = Fold[
(Print[{#[[1]], ByteCount@#[[2, 1]], #[[2, 2]]}]; #[[2, 1]]) &@
AbsoluteTiming @ {Map[Simplify, #1, {#2}], #2} &,
solnExp, -Range[6, 12, 2]]; // AbsoluteTiming
(* timing bytecount level
{18.190976, 465295824, -6}
{11.706499, 261725136, -8}
{3.647407, 225441152, -10}
{17.897693, 166903336, -12}
{52.811055, Null}
*)
Depth[soln1]
(*
36
*)
We can try to go further up, levels -15, -20, -25, -30, but there are diminishing returns:
soln2 = Fold[
(Print[{#[[1]], ByteCount@#[[2, 1]], #[[2, 2]]}]; #[[2, 1]]) &@
AbsoluteTiming@{Map[Simplify, #1, {#2}], #2} &,
soln1, -Range[15, 30, 5]]; // AbsoluteTiming
(*
{338.471629, 146285856, -15}
hmm…if it ever finishes I'll update.
*)
