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

  • $\begingroup$ Where is your try? At least your should type the equations in the mathemtica syntax. $\endgroup$
    – zhk
    Mar 25, 2016 at 4:33
  • $\begingroup$ As to the "become True" part, read this: mathematica.stackexchange.com/a/46239/1871 $\endgroup$
    – xzczd
    Mar 25, 2016 at 5:01

1 Answer 1


Even if we provide random initial conditions the DSolve is still not giving any output. So, I tried NDSolve for the numeric solution. For this I have used random values for c1, c2 and IC's

system = {l0'[t] == -l1[t], 
   l1'[t] == -l1[t] + l2[t] - c1*c2*l1[t]/(l1[t] + h1[t]), 
   l2'[t] == -l2[t] + c1*c2*l1[t]/(l1[t] + h1[t]), 
   h0'[t] == c1*(1 - c2) - h1[t];
   h1'[t] == 
    c1*(1 - c2) - h1[t] + h2[t] - c1*c2*h1[t]/(l1[t] + h1[t]), 
   h2'[t] == -h2[t] + c1*c2*h1[t]/(l1[t] + h1[t]), h1[0] == 1, 
   h2[0] == 2, l0[0] == 3, l1[0] == 4, l2[0] == 5};

c1 = 1; c2 = 1;

sol = NDSolve[system, {h1[t], h2[t], l0[t], l1[t], l2[t]}, {t, 0, 15}]

Plot[{h1[t] /. sol, h2[t] /. sol}, {t, 0, 15}, PlotRange -> All, 
 PlotStyle -> {Blue, Red}, 
 PlotLegends -> Placed[{"h1[t]", "h2[t]"}, Above]]

enter image description here

  • $\begingroup$ So if mathematica did not give any output, does this mean that mathematica cannot solve it? or the user input something wrong? or other things? For matlab, if it cannot find analytical solution, then it will generate a warning message stating explicitly that it cannot find an analytical solution. $\endgroup$
    – KevinKim
    Mar 25, 2016 at 5:38
  • $\begingroup$ I think DSolve is unable to find exact solution to this coupled non-linear system of ODEs. $\endgroup$
    – zhk
    Mar 25, 2016 at 5:41
  • $\begingroup$ I am try to solve a simpler but similar version of my problem. I type " system = {x'[t] == y[t] - c1*c2*x[t]/(1 - c1*(1 - exp (-t))), y'[t] == -y[t] + c1*c2*x[t]/(1 - c1*(1 - exp (-t)))}; " and try to solve is "sol = DSolve[system, {x, y}, t]", still Mathematica didn't produce anything. It is quite annoying, since it doesn't say that "no closed-form solution found", it just doesn't output anything. It does make my input looks nicer. But how we suppose to know that if it is because Mathematica cannot solve the problem or it's because we input something wrong? $\endgroup$
    – KevinKim
    Mar 25, 2016 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.