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

Question

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;