# A system of linear differential equations with variable (time-dependent) coefficients

I have searched almost all similar pages discussing a system of differential equations, but none of them discussed the case of variable coefficients. My problem is a bit more complicated, because instead of functions we have quantum mechanical operators (say, $$o_1(t)$$ or $$o_1^\dagger(t)$$, which is the hermitian conjugate of $$o_1(t)$$. I have shown $$\dagger$$ and $$*$$ signs by subscript $$s$$ in the following codes). One cannot assign an initial value to $$o_1$$, but it is meaningful to have expressions like $$\langle o_1^\dagger(0) o_1(0)\rangle$$, which depends on the quantum mechanical state. In other words, $$\langle o_1^\dagger(0) o_1(0)\rangle$$ is known, but $$o_1(0)$$ is unknown.

Now consider a system of differential equaitons like where

 m11 = -0.25;
m12 = 0;
m13 = 0.75 E^((0. + 22.8 I) t);
m14 = 0;
m15 = (0.097 + 0.193 E^((0. - 22.8 I) t) -
0.29 E^((0. - 45.6 I) t));
m16 = (-0.29 + 0.193 E^((0. + 22.8 I) t) +
0.097 E^((0. + 45.6 I) t));m21 = 0;
m22 = -0.25;
m23 = 0;
m24 =  0.75 E^((0. - 22.8 I) t) ; m25 = (-0.29 +
0.193  E^((0. - 22.8 I) t) + 0.0967 E^((0. - 45.6 I) t));
m26 = (0.097 + 0.193 E^((0. + 22.8 I) t) -
0.29 E^((0. + 45.6 I) t));
m31 = 0.75 E^((0. - 22.8 I) t);
m32 = 0;
m33 = -0.25;
m34 = 0;
m35 = (0.097 E^((0. - 22.8 I) t) + 0.193 E^((0. - 45.6 I) t) -
0.29 E^((0. - 68.4 I) t));
m36 = (0.193 - 0.29 E^((0. - 22.8 I) t) +
0.097 E^((0. + 22.8 I) t));
m41 = 0.75 E^((0. + 22.8 I) t) ;
m42 = 0;
m43 = -0.25;
m44 = 0;
m45 = (0.193 + 0.097 E^((0. - 22.8 I) t) -
0.29 E^((0. + 22.8 I) t));
m46 = (0.097 E^((0. + 22.8 I) t) + 0.193 E^((0. + 45.6 I) t) -
0.29 E^((0. + 68.4 I) t));
m51 = (-0.097 - 0.193 E^((0. + 22.8 I) t) +
0.29 E^((0. + 45.6 I) t));
m52 = (-0.29 + 0.193 E^((0. + 22.8 I) t) +
0.097 E^((0. + 45.6 I) t));
m53 = (-0.097 E^((0. + 22.8 I) t) - 0.193 E^((0. + 45.6 I) t) +
0.29 E^((0. + 68.4 I) t));
m54 = (0.193 - 0.29 E^((0. - 22.8 I) t) +
0.097 E^((0. + 22.8 I) t));
m55 = -0.037;
m56 = 0;
m61 = (-0.29 + 0.193 E^((0. - 22.8 I) t) +
0.097 E^((0. - 45.6 I) t));
m62 = (0.097 + 0.193 E^((0. - 22.8 I) t) -
0.29 E^((0. - 45.6 I) t)) ;
m63 = (0.193 + 0.097 E^((0. - 22.8 I) t) -
0.29 E^((0. + 22.8 I) t));
m64 = (0.097 E^((0. - 22.8 I) t) + 0.193 E^((0. - 45.6 I) t) -
0.29 E^((0. - 68.4 I) t));
m65 = 0;
m66 = -0.037;

\[Lambda]1 = (0. -
2907.509 I) + (0. +
5330.436 I) E^((0. - 22.8 I) t) - (0. +
2422.925 I) E^((0. + 22.8 I) t);

\[Lambda]1s = (0. +
2907.511 I) - (0. + 5330.437 I) E^((0. + 22.8 I) t) + (0. +
2422.925 I) E^((0. - 22.8 I) t);

\[Lambda]2 = (0. -
1453.754 I) - (0. + 4845.851 I) E^((0. - 22.8 I) t) + (0. +
6299.608 I) E^((0. - 45.6 I) t);

\[Lambda]2s = (0. +
1453.755 I) + (0. + 1356.17 I) E^((0. + 22.8 I) t) - (0. +
6299.608 I) E^((0. + 45.6 I) t);

\[Lambda]3 = (0. -
1502.866 I) - (0. +
1127.149 I) E^((0. - 22.8 I) t) + (0. +
5260.031 I) E^((0. + 22.8 I) t) - (0. +
1502.866 I) E^((0. + 45.6 I) t) - (0. +
1127.149 I) E^((0. + 68.4 I) t);

\[Lambda]3s = (0. +
1502.866 I) - (0. +
5260.031 I) E^((0. - 22.8 I) t) + (0. +
1127.149 I) E^((0. + 22.8 I) t) + (0. +
1502.866 I) E^((0. - 45.6 I) t) + (0. +
1127.149 I) E^((0. - 68.4 I) t);

n1 = +0.707 E^((0. + 34.2 I) t) fas + E^((0. + 34.2 I) t) fp;
n1s = 0.707 E^((0. - 34.2 I) t) fa + E^((0. - 34.2 I) t) fps;
n2 = 0.707 E^((0. + 11.4 I) t) fas - E^((0. + 11.4 I) t) fp;
n2s = 0.707 E^((0. - 11.4 I) t) fa - E^((0. - 11.4 I) t) fps;
n3 = 0.272 E^((0. + 22.8 I) t) fm;
n3s = 0.272 E^((0. - 22.8 I) t) fms;
M = {{m11, m12, m13, m14, m15, m16}, {m21, m22, m23, m24, m25,
m26}, {m31, m32, m33, m34, m35, m36}, {m41, m42, m43, m44, m45,
m46}, {m51, m52, m53, m54, m55, m56}, {m61, m62, m63, m64, m65,
m66}};

\[Lambda] = {{\[Lambda]1}, {\[Lambda]1s}, {\[Lambda]2}, \
{\[Lambda]2s}, {\[Lambda]3}, {\[Lambda]3s}};

n = {{n1}, {n1s}, {n2}, {n2s}, {n3}, {n3s}};


In page 8 of https://arxiv.org/pdf/1609.00075.pdf it is suggested to solve the system by the matrix exponential:

In general, it takes very long for Mathematica to calculate the matrix exponential. So, it is suggested to expand $$\mathcal T \int_{0}^{t} d\tau M(\tau)$$ as follows:

where the time interval is divided into $$N$$ pieces: $$h=(t-0)/N$$. To ensure the time ordering, $$\mathcal T$$, the lower limit of the product is $$N-1$$ and the upper is $$0$$. The reference has defined As you see, the lower limit of integral is $$\tau$$, not $$0$$. The whole problem is summarized in finding the above matrix and then calculating expressions like $$A=\int_{0}^{t}d\tau d_{11}(t,\tau) \lambda_1(\tau)$$. Finally, one can plot $$A^*A$$ versus time. I have no idea how to divide the exponential integral where the lower limit $$\tau$$ is another integration variable. As you see, we don't need to know the initial value of $$o_1$$ and so on. As one can find $$d_{ij}(t,\tau)$$ and then integrate with respect to $$\tau$$, everything is done.

Finally, I should emphasize that one can find the exponential integral analytically, but in that case one have to use MatrixExp` which for some parameters takes very long.

• So... where is the question? – Henrik Schumacher Jan 13 at 6:52
• @HenrikSchumacher I have no idea how to write a code that calculates the above equaion (exponential integral), knowing M. I have tried MatrixExp, but in general it doesn't work well. – Saeid Jan 13 at 14:06
• Your code is lacking the first row of $M$. – Cesareo Jan 13 at 18:17
• @Cesareo I fixed that. – Saeid Jan 13 at 20:29