I have a set of coupled non-linear ODEs, they look symmetric. But when I use matlab symbolic calculation to solve, matlab told me that no explicit solution found. I was recommended by my friend to try mathematica. But I'm not very familiar with the syntax in mathematica. I followed the example in online document of mathematica but it seems does not work for my case.
The set of ODE is the following (let $q\in(0,1)$ and $\rho\in(0,1)$ be constants): \begin{equation} \begin{cases} &\dot{\pi}^L_0(t) = \pi^L_1(t)\\ &\dot{\pi}^L_1(t) = -\pi^L_1(t)+\pi^L_2(t)-\rho(1-q)\frac{\pi^L_1(t)}{\pi^L_1(t)+\pi^H_1(t)}\\ &\dot{\pi}^L_2(t) = -\pi^L_2(t)+\rho(1-q)\frac{\pi^L_1(t)}{\pi^L_1(t)+\pi^H_1(t)} \end{cases} ~\text{and }~\begin{cases} &\dot{\pi}^H_0(t) = -\rho q+\pi^H_1(t)\\ &\dot{\pi}^H_1(t) = \rho q-\pi^H_1(t)+\pi^H_2(t)-\rho(1-q)\frac{\pi^H_1(t)}{\pi^L_1(t)+\pi^H_1(t)}\\ &\dot{\pi}^H_2(t) = -\pi^H_2(t)+\rho(1-q)\frac{\pi^H_1(t)}{\pi^L_1(t)+\pi^H_1(t)} \end{cases} \end{equation} The initial condition is $\pi^L_1(0)>0$, $\pi^L_0(0)>0$ and $\pi^L_2(0)=0$ with $\pi^L_0(0)+\pi^L_1(0)+\pi^L_2(0)=p_L\in(0,1)$ and $\pi^H_1(0)>0$, $\pi^H_0(0)>0$ and $\pi^H_2(0)=0$ with $\pi^H_0(0)+\pi^H_1(0)+\pi^H_2(0)=1 - p_L$. Note that I can define $\pi_j(t)=\pi^L_j(t)+\pi^H_j(t)$, then I can solve $\pi_j(t)$ analytically, which is \begin{equation} \begin{cases} &\pi_2(t)=\rho(1-q)(1-e^{-t})\\ &\pi_1(t)=\rho q-\rho(1-q)e^{-t}t+(\pi_1(0)-\rho q)e^{-t}\\ &\pi_0(t)=1-\pi_1(t)-\pi_2(t) \end{cases} \end{equation} And observe that $\dot{\pi}^L_0+\dot{\pi}^L_1+\dot{\pi}^L_2=0$ and $\dot{\pi}^H_0+\dot{\pi}^H_1+\dot{\pi}^H_2=0$, this implies that $\pi^L_0(t)+\pi^L_1(t)+\pi^L_2(t)=p_L$ and $\pi^H_0(t)+\pi^H_1(t)+\pi^H_2(t)=p_L$.
Note that after observing these fact, the H-part ODE and L-part ODE can be decoupled. So we only need to solve H-part and L-part separately. However, even if take these fact into consideration, Mathematica still cannot solve it I typed the following command
system = {x'[t] == y[t],
y'[t] == -y[t] + z[t] - c1*c2*y[t]/(c1*(1 - c2) -
c1*c2*Exp[-t]*t + (c3 - c1*(1 - c2))*Exp[-t]),
z'[t] == -z[t] + c1*c2*y[t]/(c1*(1 - c2) -
c1*c2*Exp[-t]*t + (c3 - c1*(1 - c2))*Exp[-t])};
here $x=\pi^L_0, y=\pi^L_1,z=\pi^L_2$
Then I type
DSolve[system, {x[t], y[t], z[t]}, t]
Mathematica only makes my input looks nicer but did not solve this problem at all
