About solving a system of first order ODEs with non constant coefficients

Question

I'm trying to solve this system of ODEs \begin{align*} &S'(t)=-\lambda(t)S(t)+{\tau}R(t)\\ &I'(t)=\lambda(t)S(t)-{\gamma}I(t)+r(t)L(t)\\ &R'(t)=(1-q){\gamma}I(t)-{\tau}R(t)\\ &L'(t)=q{\gamma}I(t)-r(t)L(t) \\ \end{align*} so if I write it in compact form: $X'(t)=A(t)X(t)$. Now, my problem is the matrix A, because it also depends on t. I have read that in order to solve this problem I should use exponential matrix, but in order to solve it I have to prove that the matrix commute with its integral i.e. $$ A(t)\int_{0}^{t}A(s)ds = \int_{0}^{t}A(s)ds{A(t)} $$ But unfortunately it does not. So is there another procedure that I can try to use to solve this system analitically or is it impossible and I am force to find a numerical one? I hope that someone of you can give me some hints-suggestions. Thank you in advance!