I have a system of differential equations that I have reduced to this:

\begin{equation} \small \frac{db}{dt}= \beta_1 (1-b(t)-c(t)) b(t) + (1-\beta_1) \beta_3 (1-b(t)-c(t)) d(t) - \phi_{w1} b(t) c(t) - \phi_{b1} b(t) e(t) - \mu_B b(t) - \gamma_1 b(t) \end{equation} \begin{equation} \small \frac{dc}{dt}= \beta_2 (1-b(t)-c(t)) c(t) + (1-\beta_2) \beta_4 (1-b(t)-c(t)) e(t) + \phi_{w1} b(t) c(t) + \phi_{b1} b(t) e(t) - \mu_C c(t) - \gamma_2 c(t) \end{equation} \begin{equation} \small \frac{dd}{dt}= \beta_3 (1-d(t)-e(t)) d(t) + (1-\beta_3) \beta_1 (1-d(t)-e(t)) b(t) - \phi_{w2} d(t) e(t) - \phi_{b2} d(t) c(t) - \mu_D d(t) - \gamma_3 d(t) \end{equation} \begin{equation} \small \frac{de}{dt}= \beta_4 (1-d(t)-e(t)) e(t) + (1-\beta_4) \beta_2 (1-d(t)-e(t)) c(t) + \phi_{w2} d(t) e(t) + \phi_{b2} d(t) c(t) - \mu_E e(t) - \gamma_4 e(t) \end{equation}

Here the code for the equations as input in mathematica:

    ode1 = b'[t] == 
   Subscript[\[Beta], 1]*(1 \[Minus] b[t] \[Minus] c[t])*
     b[t] + (1 \[Minus] Subscript[\[Beta], 1])*
     Subscript[\[Beta], 3]*(1 \[Minus] b[t] \[Minus] c[t])*
     d[t] \[Minus] Subscript[\[Phi], w1]*b[t]*c[t] \[Minus] 
    Subscript[\[Phi], b1]*b[t]*e[t] \[Minus] 
    Subscript[\[Mu], B]*b[t] \[Minus] Subscript[\[Gamma], 1]*b[t];
ode2 = c'[t] == 
   Subscript[\[Beta], 2]*(1 \[Minus] b[t] \[Minus] c[t])*
     c[t] + (1 \[Minus] Subscript[\[Beta], 2])*
     Subscript[\[Beta], 4]*(1 \[Minus] b[t] \[Minus] c[t])*e[t] + 
    Subscript[\[Phi], w1]*b[t]*c[t] + 
    Subscript[\[Phi], b1]*b[t]*e[t] \[Minus] 
    Subscript[\[Mu], C]*c[t] \[Minus] Subscript[\[Gamma], 2]*c[t];
ode3 = d'[t] == 
   Subscript[\[Beta], 3] *(1 \[Minus] d[t] \[Minus] e[t])*
     d[t] + (1 \[Minus] Subscript[\[Beta], 3])*
     Subscript[\[Beta], 1]*(1 \[Minus] d[t] \[Minus] e[t])*
     b[t] \[Minus] Subscript[\[Phi], w2]*d[t]*e[t] \[Minus] 
    Subscript[\[Phi], b2]*d[t]*c[t] \[Minus] 
    Subscript[\[Mu], D]*d[t] \[Minus] Subscript[\[Gamma], 3]*d[t];
ode4 = e'[t] == 
   Subscript[\[Beta], 4] *(1 \[Minus] d[t] \[Minus] e[t])*
     e[t] + (1 \[Minus] Subscript[\[Beta], 4])*
     Subscript[\[Beta], 2]*(1 \[Minus] d[t] \[Minus] e[t])*c[t] + 
    Subscript[\[Phi], w2]*d[t]*e[t] + 
    Subscript[\[Phi], b2]*d[t]*c[t] \[Minus] 
    Subscript[\[Mu], E]*e[t] \[Minus] Subscript[\[Gamma], 4]*e[t];

Would the code to solve the system in mathematica be something like this?

  b[0] = 100
  c[0] = 500
  d[0] = 100
  e[0] = 400
  (fb, fc, fd, fe} = 
 NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == k1, c[0] == k2, 
    d[0] == k3, e[0] == k4}, {b, c, d, e}, {t, 0, 1000}].

I was trying to generate a visualisation of the Phase portrait of this system in mathematica.

Here b, c, d and e are positive real variables (population). Variable t goes from 0 to 1000. All the rest of parameters are all decimal numbers taking values from 0 to 1. We can assume constant prameters $\mu$ = 0.06 and constant parameters $\gamma$=0.01 in all the cases.

How can I a make a good visualisation of phase portraits of this system. For example, a 3D plot.

Can anyone help me?

  • 8
    $\begingroup$ At a minimum, edit your question to include the Mathematica code for your equations in copy and paste-able form (InputForm). Make all dependencies on t explicit. $\endgroup$
    – Bob Hanlon
    Commented Mar 25, 2022 at 16:04
  • 3
    $\begingroup$ Aside from the comment by @BobHanlon which is important, I'd like to add that you present 4 equations, hence, naively, we are led to believe you want a 4-dimensional plot which -as far as I know- is not possible. After presenting your equations, you mention that you want a 3 dimensional visualization, which means that a projection to 3 dimensions is needed. This makes sense because the dimensions of a phase space is equal to how many independent variables you have. Question: what are the independent variables? Could you re-write the system in a better way? $\endgroup$
    – bmf
    Commented Mar 25, 2022 at 18:43
  • 3
    $\begingroup$ @pyring: Would be a lot easier to read if you just used $y,z,w,u$ for dependent variables and $a,b,c,d $ and so forth for constants (the dd notation is unnecessarily confussing). Also a lot easier to input into Mathematica that way. Maybe you should review `NDSolve' help (type it in , hover over it, press "I", and study some examples. Then try and input your system with simple initial conditions to get NDSolve working and producing a result. Then paste your code back here to go further even if it's not working but can begin there. $\endgroup$
    – josh
    Commented Mar 25, 2022 at 19:54
  • 1
    $\begingroup$ You need a * or a space between variables that you multiple (e.g. bc) otherwise Mathematica assumes that bc is a new variable (not b*c). $\endgroup$
    – Chris K
    Commented Mar 26, 2022 at 19:45
  • 2
    $\begingroup$ You need to establish initial conditions b(0)=b0,c(0)=c0,d(0)=d0,e(0)=e0 and you need to set all the constants to some initially-arbitrary values just to get NDSolve (or NDSolveValue) up and running without errors. Make these and the changes above, and try and set up your system to run: (fb,fc,fd,fe}=NDSolveValue[{ode1,ode2,ode3,ode4,b[0]==k1,c[0]==k2,d[0]==k3,e[0]==k4},{b,c,d,e},{t,0,1000}]. Post your code. $\endgroup$
    – josh
    Commented Mar 26, 2022 at 20:07

1 Answer 1


It's a good try pyring. But in order to solve first-order ODEs numerically, you need two things: (1) Supply numeric values to all constants, (2) Supply numeric values for initial conditions. You have 16 constants and four equations so that's 20 values you need to assign. In the code below I chose random values between 0 and 1. The important thing is to just get it running without syntax errors or run-time errors. Later you can adjust it to your needs. In the code below I plotted individually the four functions written as fb[t], fc[t],fd[t] and fe[t]. Then I chose three of them fb[t], fc[t] and fd[t] to construct a 3D phase portrait which appears as a simple line. Maybe with other choices of initial conditions and constants you can achieve something more interesting. Just cut this code and paste it into your notebook and it should run without errors. Then you can experiment with changing the initial conditions and constants.

     set up constants and ODEs
    Subscript[\[Mu], B] = 0.5;
     Subscript[\[Beta], 1] = 0.1;
    Subscript[\[Beta], 3] = 0.2;
    Subscript[\[Gamma], 1] = 0.7;
    Subscript[\[Gamma], 2] = 0.221;
    Subscript[\[Phi], b1] = 0.9;
    Subscript[\[Phi], w1] = 0.05;
    Subscript[\[Phi], w2] = 0.08;
    Subscript[\[Beta], 2] = 0.8;
    Subscript[\[Beta], 4] = 0.2;
    Subscript[\[Mu], C] = 0.5;
     Subscript[\[Mu], E] = 0.01;
    Subscript[\[Phi], b2] = 0.6;
    Subscript[\[Gamma], 3] = 0.109;
    Subscript[\[Mu], D] = .07;
    Subscript[\[Gamma], 4] = 0.219;
    ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) + 
        b[t] Subscript[\[Beta], 
         1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 
           1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) - 
        b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] - 
        b[t] c[t] Subscript[\[Phi], w1];
    ode2 = Derivative[1][c][t] == 
       Subscript[\[Beta], 2]c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta],2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) + 
        b[t] e[t] Subscript[\[Phi], b1] + 
        b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] - 
        c[t] Subscript[\[Mu], C];
    ode3 = Derivative[1][d][t] == 
       b[t] Subscript[\[Beta], 
         1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) - 
        d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] - 
        d[t] Subscript[\[Mu], D] + 
        Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) - 
        d[t] e[t] Subscript[\[Phi], w2];
    ode4 = Derivative[1][e][t] == 
       d[t] c[t] Subscript[\[Phi], b2] + 
         2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) + 
        Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) + 
        d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] - 
        e[t] Subscript[\[Mu], E];
     run system with some arbitrary initial conditions
    {fb, fc, fd, fe} = 
      NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2, 
        d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}];
     plot each function separately
    Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000}, 
     PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]]
     choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \
    ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100}, 
     BoxRatios -> {1, 1, 1}, 
     PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]

enter image description here

enter image description here

  • $\begingroup$ thanks. However, I copy paste the code but I obtain a different plot in Mathematica 11.3. $\endgroup$
    – pyring
    Commented Mar 28, 2022 at 15:42
  • 1
    $\begingroup$ Ok. Not sure why. I just tried it again, copy, paste into a new notebook and it gives the same result. You can try separating the odes into separate cells and make sure they are formatted correctly: Highlight just ode1 for example, then press Shift/Ctrl D. This will create a separate cell for just theat ode. Now press Enter and see if the ode is as you expected it. Then do the same with others. Then run NDSolveValue again. $\endgroup$
    – josh
    Commented Mar 29, 2022 at 11:33

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.