I try to solve a system of 16 linear differential equations, with a following code:
sminus = {{0, 0}, {1, 0}};
splus = {{0, 1}, {0, 0}};
kappa[t_, n_, d1_,g_] = {{xk[t, n, d1, g], yk[t, n, d1, g]}, {yk1[t, n, d1, g],zk[t, n, d1, g]}};
kappam[t_, n_, d1_,g_] = {{xkm[t, n, d1, g], ykm[t, n, d1, g]}, {ykm1[t, n, d1, g],zkm[t, n, d1, g]}};
kappamp[t_, n_, d1_,g_] = {{xkmp[t, n, d1, g], ykmp[t, n, d1, g]}, {ykmp1[t, n, d1, g],zkmp[t, n, d1, g]}};
kappapm[t_, n_, d1_,g_] = {{xkpm[t, n, d1, g], ykpm[t, n, d1, g]}, {ykpm1[t, n, d1, g],zkpm[t, n, d1, g]}}
DSolve[Flatten@{D[kappa[t, n, d1, g], t] ==
-n^2*(sminus.splus.kappa[t, n, d1, g] +
kappa[t, n, d1, g].sminus.splus) - (n +
1)^2*(splus.sminus.kappa[t, n, d1, g] +
kappa[t, n, d1, g].splus.sminus) +
2*n*(n + 1)*(splus.kappa[t, n, d1, g].sminus +
sminus.kappa[t, n, d1, g].splus
-
2*n*(ConjugateTranspose[kappam[t, n, d1, g]].sminus +
splus.kappam[t, n, d1, g]) +
2*(n + 1)*(kappam[t, n, d1, g].splus +
sminus.ConjugateTranspose[kappam[t, n, d1, g]]) +
2*n*kappamp[t, n, d1, g] - 2*(n + 1)*kappapm[t, n, d1, g]) +
I*(sminus.splus.kappa[t, n, d1, g] -
kappa[t, n, d1, g].sminus.splus)*d1,
D[kappam[t, n, d1, g], t] ==
-n^2*(sminus.splus.kappam[t, n, d1, g] +
kappam[t, n, d1, g].sminus.splus) - (n +
1)^2*(splus.sminus.kappam[t, n, d1, g] +
kappam[t, n, d1, g].splus.sminus) +
2*n*(n + 1)*(splus.kappam[t, n, d1, g].sminus +
sminus.kappam[t, n, d1, g].splus)
- n^2*kappam[t, n, d1, g] - (n + 1)^2*kappam[t, n, d1, g] +
2*((n + 1)^2*sminus.kappapm[t, n, d1, g] +
n^2*kappamp[t, n, d1, g].sminus -
n*(n + 1)*sminus.kappamp[t, n, d1, g] -
n*(n + 1)*kappapm[t, n, d1, g].sminus)
+ I*(sminus.splus.kappam[t, n, d1, g] -
kappam[t, n, d1, g].sminus.splus)*d1 -
I*d1*kappam[t, n, d1, g],
D[kappamp[t, n, d1, g], t] ==
-n^2*(sminus.splus.kappamp[t, n, d1, g] +
kappamp[t, n, d1, g].sminus.splus) - (n +
1)^2*(splus.sminus.kappamp[t, n, d1, g] +
kappamp[t, n, d1, g].splus.sminus) +
2*n*(n + 1)*(splus.kappamp[t, n, d1, g].sminus +
sminus.kappamp[t, n, d1, g].splus)
+ 2*n*(n + 1)*kappapm[t, n, d1, g] -
2*n^2*kappamp[t, n, d1, g] +
2*((n + 1)^2*
sminus.ConjugateTranspose[
kappam[t, n, d1, g]] + (n + 1)^2*
kappam[t, n, d1, g].splus -
n*(n + 1)*ConjugateTranspose[kappam[t, n, d1, g]].sminus -
n*(n + 1)*splus.kappam[t, n, d1, g])
+ I*(sminus.splus.kappamp[t, n, d1, g] -
kappamp[t, n, d1, g].sminus.splus)*d1,
D[kappapm[t, n, d1, g], t] ==
-n^2*(sminus.splus.kappapm[t, n, d1, g] +
kappapm[t, n, d1, g].sminus.splus) - (n +
1)^2*(splus.sminus.kappapm[t, n, d1, g] +
kappapm[t, n, d1, g].splus.sminus) +
2*n*(n + 1)*(splus.kappapm[t, n, d1, g].sminus +
sminus.kappapm[t, n, d1, g].splus)
+ 2*n*(n + 1)*kappamp[t, n, d1, g] -
2*(n + 1)^2*kappapm[t, n, d1, g] +
2*(n^2*ConjugateTranspose[kappam[t, n, d1, g]].sminus +
n^2*splus.kappam[t, n, d1, g] -
n*(n + 1)*sminus.ConjugateTranspose[kappam[t, n, d1, g]] -
n*(n + 1)*kappam[t, n, d1, g].splus)
+ I*(sminus.splus.kappapm[t, n, d1, g] -
kappapm[t, n, d1, g].sminus.splus)*d1},
Flatten@{kappa[t, n, d1, g], kappam[t, n, d1, g],
kappamp[t, n, d1, g], kappapm[t, n, d1, g]}, t]
Unfortunately Mathematica just instantly returns back the input without even attempting to solve it. It is not so relevant for me how long it would take to calculate the solution, but I would rather like to avoid to use NDSolve if it is possible.