# Help speeding up matrix exponent of $16 \times 16$ symbolic matrix

I am trying to run the following script which involves the MatrixExp of t times a 16x16 matrix with one variable, v, in order to solve a 16x16 system of coupled differential equations.

h2 = {{0, 0, -Ω1, 0}, {0, -Δ1 + Δ2 + k1 v + k2 v, -Ω2, 0}, {-Ω1, -Ω2, -Δ1 +k1 v, 0}, {0, 0, 0, 0}};
ρ = {{ρ11[t], ρ12[t], ρ13[t], 0}, {ρ21[t], ρ22[t], ρ23[t], 0}, {ρ31[t], ρ32[t], ρ33[t], 0}, {0, 0, 0, 0}};
ρvar = {{ρ11[t], ρ12[t], ρ13[t], ρ14[t]}, {ρ21[t], ρ22[t], ρ23[t], ρ24[t]}, {ρ31[t], ρ32[t], ρ33[t], ρ34[t]}, {ρ41[t], ρ42[t], ρ43[t], ρ44[t]}};
ρprime = -I (h2.ρ - ρ.h2) + {{γ31 ρ33[t] + γ21 ρ22[t], -(1/2) γ21 ρ12[t], -(1/2) (γ31 + γ32) ρ13[t], 0}, {-(1/2) γ21 ρ21[t], -γ21 ρ22[t] + γ32 ρ33[t], -(1/2) (γ21 + γ31 + γ32) ρ23[t], 0}, {-(1/2) (γ31 + γ32) ρ31[t], -(1/2) (γ21 + γ31 + γ32) ρ32[t], -ρ33[t] (γ31 + γ32 + γ34), 0}, {0, 0, 0, ρ33[t] γ34}};
replace3 = {Δ1 -> (2.1 π)/(500*10^-7)*10^3,
Δ2 -> (2 π)/(500*10^-7)*10^3, γ21 ->
1/(16*10^-9), γ31 -> 1/(16*10^-9), γ32 ->
1/(16*10^-9), γ34 -> 0.1/(16*10^-9), Ω1 ->
10^9, Ω2 -> 10^9, k1 -> (2.1 π)/(500*10^-7),
k2 -> (2 π)/(500*10^-7), m -> 10^-25, ℏ -> 10^-34,
c -> 3*10^8};

var = Flatten@ρvar;

{barray, marray} =
CoefficientArrays[Flatten@ρprime /. replace3, var];
presol = MatrixExp[
marray t, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}] //
AbsoluteTiming;
sol = Flatten@presol;

• You could try to compute a symbolic JordanDecomposition first. In in general, working with symbolic matrix for size greater than $4 \times 4$ should be avoided (unless they have a special structure that can be exploited). – Henrik Schumacher Feb 25 at 8:44
• (1) Could try making the input exact since the mix of approximate numbers and symbolic values is problematic to symbolic computation. (2) Alternatively, define a function that only evaluates for numeric values of v. (3) All that aside, it might be slow due to a problem in Eigenvalues that I will investigate. (Those suggestions will bypass the problem, whether they turn out to be useful or not.) – Daniel Lichtblau Feb 25 at 17:47
• I recommend you write down the actual differential equations and then use NDSolve to solve them. Alternatively, you could look into the Magnus expansion to get good approximations to the analytic solution. – Roman Feb 26 at 17:53