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.

*)