plot system of 5 odes over time using ode45

8 views (last 30 days)
Sam Rofman on 8 Dec 2021
Edited: Walter Roberson on 11 Dec 2021
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;

Walter Roberson on 8 Dec 2021
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 on 11 Dec 2021
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