plot system of 5 odes over time using ode45

Here are my equations and initial conditions. I'm trying to plot these odes as a function of time with the following initial conditions:

xinit = [10000, 0, 0.01, 0, 0];
f = @(t,y) [(10000 - (0.01*y(1)) - (0.000000024*y(3)*y(1)*20001)); ((0.05*0.000000024*y(3)*y(1)*20001) - y(2)); ((25000*y(2)) - 23); ((0.95*0.000000024*y(3)*y(1)*20001) - (0.001*y(4))); ((15*0.001*y(4)) - (6.6*y(5)))];

I'm looking for oscillatory behavior. Here they are written slightly differently:

      T = 10000;
      Q = 0;
      V = 0.01
      P = 0;
      I = 0;
      dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;
      dydt(2) = 0.05*0.000000024*V*T*20001 - Q;
      dydt(3) = 25000*Q - 23;
      dydt(4) = 0.95*0.000000024*V*T*20001 - 0.001*P;
      dydt(5) = 15*0.001*P - 6.6*I;

Answer 1

[Walter Roberson]

tspan = [0 0.184];
xinit = [10000, 0, 0.01, 0, 0];
[t, y] = ode45(@odefun, tspan, xinit);

plot(t, y);
legend({'T', 'Q', 'V', 'P', 'I'}, 'location', 'best');

function dydt = odefun(t, y)
    T = y(1); Q = y(2); V = y(3); P = y(4); I = y(5);
    
    dydt = zeros(5,1);

    dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;
    dydt(2) = 0.05*0.000000024*V*T*20001 - Q;
    dydt(3) = 25000*Q - 23;
    dydt(4) = 0.95*0.000000024*V*T*20001 - 0.001*P;
    dydt(5) = 15*0.001*P - 6.6*I;
end

You cannot go much further. By roughly 1.88 the ode functions refuse to continue as the slopes have become to steep.

[Walter Roberson]

Look at the first equation, and compare to

    dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;

The equation starts with the constant s, which matches the 10000 in your code; so far so good.

Then the equation has - and we can speculate that the maps to your constant named dt (not a great name considering that you are taking derivatives with respect to t and so have dt in your equations with quite a different meaning... but anyhow.) That plausibly matches the - 0.01*T term, so that's okay.

Then the equation has . k is a constant, V is the function we have mapped to y(3), T is the function we have mapped to y(1), so the part up to 0.000000024*V*T matches. But after that in the equation is the part. γ is a constant. But I is the function mapped to y(5) . So the code should not have 20001 at that point, and should instead have (1 + 2*10^4 * I) where I = y(5) is in effect.

All of your 20001 in your code should be (1 + 2*10^4 * I)

You should re-check the code.

And there is no shame at all in putting all of those constants such as s as variables in your function code and using them to make the code clear, instead of hard-coding values like 0.000000024 or 20000 in the equations. Any tiny loss you might take in execution speed for using more variables will be more than made up for in time saved debugging the code.

[Walter Roberson]

tspan = [0 1];
xinit = [10000, 0, 0.01, 0, 0];
[t, y] = ode45(@odefun, tspan, xinit);

symbols = {'T', 'T^*', 'V', 'P', 'I'};
N = size(y,2);

for K = 1 : N
  subplot(N,1,K);
  plot(t, y(:,K));
  legend(symbols(K), 'location', 'best');
end

function dydt = odefun(t, y)
    T = y(1); T__star = y(2); V = y(3); P = y(4); I = y(5);

    s = 10000;
    k = 2.4*10^-8;
    gamma  = 2*10^4;
    f = 0.95;
    r = 2.5*10^4;
    B = 15;
    d_T = 0.01;
    d_T__star = 1;
    d_V = 23;
    d_P = 6.6;
    d_I = 0.001;

    dydt = zeros(5,1);

    dydt(1) = s - d_T*T - k*V*T*(1+gamma*I);
    dydt(2) = (1-f)*k*V*T*(1+gamma*I) - T__star;
    dydt(3) = r*T__star - d_V;
    dydt(4) = f*k*V*T*(1+gamma*I) - d_P*P;
    dydt(5) = B*d_P*P - d_I*I;
end

[Walter Roberson]

You can use tiledlayout() to get larger graphs.

Or just put them into different figures.

Plots 2, 4, 5 look like they might fit well on the same plot, but plots 1 and 3 have a significantly different range and do not fit on the same plot.

Good thing you caught my mistake in the second equation.

[Walter Roberson]

tspan = [0:.2:50, 100:100:2000];
xinit = [1000, 0, 0.001, 0, 0]; %OR ;xinit = [1000, 0, 0.001, 0, 0]; OR ;xinit = [1, 0, 0.001, 0, 0]; OR ;xinit = [10000, 0, 0.001, 0, 0];
[t, y] = ode45(@odefun, tspan, xinit);
symbols = {'T', 'T^*', 'V', 'P', 'I'};

figure()
whichplot = [1 3];
N = length(whichplot);

for K = 1: N
  subplot(N,1,K);
  plot(t, y(:,whichplot(K)));
  legend(symbols(K), 'location', 'best');
end

function dydt = odefun(t, y)
    T = y(1); T__star = y(2); V = y(3); P = y(4); I = y(5);
    s = 10000;
    k = 2.4*10^-8;
    gamma  = 2*10^-4;
    f = 0.95;
    r = 2.5*10^4;
    B = 15;
    d_T = 0.01;
    d_T__star = 1;
    d_V = 0.001;
    d_P = 23;
    d_I = 6.6;
    dydt = zeros(5,1);
    dydt(1) = s - d_T*T - k*V*T*(1+gamma*I);
    dydt(2) = (1-f)*k*V*T*(1+gamma*I) - d_T__star*T__star;
    dydt(3) = r*T__star - d_V;
    dydt(4) = f*k*V*T*(1+gamma*I) - d_P*P;
    dydt(5) = B*d_P*P - d_I*I;
end