How to plot solutions of system of nonlinear differential equations in MATLAB?

Question

I am trying to follow this research paper. I am trying to duplicate the graph of the solutions located in Figure 7 on page 20. I have a screenshot of Figure 7:

I would first like to recreate the left picture. The system in question is what I have as dX. Here is what I have in an m-file:

function dX = CompetitionModel(t,X)       
  bs = 8*10^(-3); 
  bl = 4*10^(-3); 
  bh = 6.4*10^(-3);
  N = bs + bl + bh;
  K = 10^8;
  m1 = 2*10^(-5); 
  m2 = 9*10^(-9); 
  p = 5*10^(-13); 
  I = 10^(-3); 
  T = 0; 
  a = 0; 
  dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3));
       X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3));
       X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))];
end

ode45 has the syntax: [T,Y] = solver(odefun,tspan,y0). I get the tspan from the picture I posted. My initial conditions are:S0 = 10^4; Rl0 = 0; Rh0 = 0, so this is what I have as y0. I type the following into the command window:

>>[t,X1] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]);
>>[t,X2] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]);
>>[t,X3] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]);

MATLAB has been busy for the past 30 minutes and my laptop is starting to get really hot. So I can't get to plotting until it's finished and I don't know if there are any errors in my code. I was wondering if there was perhaps a better way that I could get the solutions of the system dX.

Answer 1

I checked the ODE against the paper and found one mistake. According to the last paragraph of page 9 in the paper, N = S + Rl + Rh. A corrected code:

function dX = CompetitionModel(~,X)       
  bs = 8e-3; 
  bl = 4e-3; 
  bh = 6.4e-3;
  N = sum(X);  % this line was incorrect
  K = 1e8;
  m1 = 2e-5; 
  m2 = 9e-9; 
  p = 5e-13; 
  I = 1e-3; 
  T = 0;
  a = 0; 
  dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3));
       X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3));
       X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))];

Note a couple of syntactic changes. First, since you do not actually use the time value in your ODE, you can leave a ~ in place of the t that you had in your function definition as a stand-in for the unused input. Second, you can use the notation 8e-3 instead of 8*10^(-3). They evaluate to the same thing, but the former looks a little cleaner.

The plot for no treatment, i.e. T = 0, is shown below.

I figured out the issue in the ODE by first trying to solve your equation with some of the stiff ODE solvers in MATLAB. See the MATLAB ODE solver documentation for more details. Basically I solved it with ode15s and ode23s and found that the solution was unstable (population went off to infinity). These other solvers are good tools to keep in mind if your solution is hanging; sometimes one of the others will work and it will either give you what you want, or indicate that you have another problem somewhere else.

Note: I am fairly certain that you are not using "phase portrait" correctly here. You are just looking for the solution of the ODE with respect to time with a given set of initial conditions. A phase portrait looks at how the states of the system (your three different populations here) evolve given different sets of initial conditions. It does not show you anything about the time-dependence of these solutions, just how they evolve relative to one-another.