How do I store an ODE in a vector?

I have a fully coded ODE modelling population growth of cells and I want to find the homeostasis levels or when the graph flattens out. in order to do this I need to store the ODEs, (tmp1, tmp2, tmp3, tmp4) which are under dydt, in vectors. How can I store each ODE into a vector and if this is confusing how is an easier way to find the point when the curve flattens and there is no change (slope=0).

Also note I cannot use the symbolic set as my license for matlab doesn't allow it.

function [t,dydt1] = HematopoieticSystemModel

time = 0:4000;

%------%

x_0 = 1;

x_1 = 3;

x_2 = 9;

x_3 = 50;

p0_0=100;

%----------%

p = [x_0 x_1 x_2 x_3 p0_0];

%--------------------------------------------------%

[t,dydt1] = ode45(@myODE,time,p);

%% -----Individual Plots------ %

figure (1)

hold on

subplot (2,2,1);

plot(t,dydt1(:, 1), 'r');

title ('HSC');

xlabel('time')

ylabel('cells')

hold off

hold on

subplot (2,2,2);

plot (t, dydt1(:,2),'g');

title ('ST-HSCs');

xlabel('time')

ylabel('cells')

hold off

hold on

subplot (2, 2, 3);

plot (t, dydt1(:,3),'b');

title ('MPPs');

xlabel('time')

ylabel('cells')

hold off

hold on

subplot (2, 2, 4);

plot (t, dydt1(:,4), 'm');

title ('CLPs and CMPs');

xlabel('time')

ylabel('cells')

hold off

%------------------%

figure (2)

hold on

plot(t,dydt1(:, 1), 'r');

plot (t, dydt1(:,2),'g');

plot (t, dydt1(:,3),'b');

plot (t, dydt1(:,4), 'm');

legend ('HSC' , 'ST-HSC', 'MPP', 'CLP and CMP');

xlabel('time (days)');

ylabel('cells');

title ('Hematopoietic System Growth');

set(gca, 'YScale', 'log');

hold off

%------------------%

figure (3)

hold on

plot(t,dydt1(:, 1), 'r');

title ('HSC');

xlabel('time')

ylabel('cells')

hold off

figure(4)

hold on

plot (t, dydt1(:,2),'g');

title ('ST-HSCs');

xlabel('time')

ylabel('cells')

hold off

figure (5)

hold on

plot (t, dydt1(:,3),'b')

title ('MPPs');

xlabel('time')

ylabel('cells')

hold off

figure (6)

hold on

plot (t, dydt1(:,4), 'm')

title ('CLPs and CMPs');

xlabel('time')

ylabel('cells')

hold off

% figure(7)

% hold on

% plot (t, dydt1(:,5), 'k')

% title ('Self Renewal for HSCs');

% xlabel('time')

% ylabel('cells')

% hold off

end

function dydt = myODE(t,y, r);

P0=0.7;

h0=0.0000235;

r0=0.01818;

r1=0.08712;

r2=8.0217;

d3=0.1;

p1=0.4783;

p2=0.4987;

%%

% y(1) = x_0;

% y(2) = x_1;

% y(3) = x_2;

% y(4) = x_3;

p0=P0./(1+h0.*y(1));

tmp1 = r0.*y(1).*(2.*p0-1);

tmp2 = 2.*r0.*y(1).*(1-p0) + r1.*y(2).*(2.*p1-1);

tmp3 = 2.*r1.*y(2).*(1-p1) + r2.*y(3).*(2.*p2-1);

tmp4 =2.*r2.*y(3).*(1-p2) - d3.*y(4);

dydt = [tmp1; tmp2; tmp3; tmp4; p0];

end

Answer 1

[Torsten]

Either you use the code below and work with the time derivatives after the integration to decide when the curves have flattened or you use the Event facility of the ODE integrators to stop integration when the time derivatives of the equations become small enough for your purpose.

%function [t,dydt1] = HematopoieticSystemModel
time = 0:4000;
%------%
x_0 = 1;
x_1 = 3;
x_2 = 9;
x_3 = 50;
p0_0=100;
%----------%
p = [x_0 x_1 x_2 x_3 p0_0];
%--------------------------------------------------%
[t,y] = ode15s(@myODE,time,p);
dydt = zeros(size(y));
for i = 1:numel(t)
  dydt(i,:) = myODE(t(i),y(i,:));
end
plot(t,dydt)
%% -----Individual Plots------ %
figure (1)
hold on

subplot (2,2,1);
plot(t,y(:, 1), 'r');
title ('HSC');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2,2,2);
plot (t, y(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 3);
plot (t, y(:,3),'b');
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 4);
plot (t, y(:,4), 'm');
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off

%------------------%
figure (2)
hold on
plot(t,y(:, 1), 'r');
plot (t, y(:,2),'g');
plot (t, y(:,3),'b');
plot (t, y(:,4), 'm');
legend ('HSC' , 'ST-HSC', 'MPP', 'CLP and CMP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off

%------------------%
figure (3)
hold on
plot(t,y(:, 1), 'r'); 
title ('HSC');
xlabel('time')
ylabel('cells')
hold off

figure(4)
hold on
plot (t, y(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off

figure (5)
hold on
plot (t, y(:,3),'b')
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off

figure (6)
hold on
plot (t, y(:,4), 'm')
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off 

% figure(7)
% hold on
% plot (t, y(:,5), 'k')
% title ('Self Renewal for HSCs');
% xlabel('time')
% ylabel('cells')
% hold off
%end
function dydt = myODE(t,y)
P0=0.7; 
h0=0.0000235;
r0=0.01818;
r1=0.08712;
r2=8.0217; 
d3=0.1;
p1=0.4783;
p2=0.4987; 
%%
% y(1) = x_0;
% y(2) = x_1;
% y(3) = x_2;
% y(4) = x_3;
p0=P0./(1+h0.*y(1));
tmp1 = r0.*y(1).*(2.*p0-1);
tmp2 = 2.*r0.*y(1).*(1-p0) + r1.*y(2).*(2.*p1-1);
tmp3 = 2.*r1.*y(2).*(1-p1) + r2.*y(3).*(2.*p2-1);
tmp4 =2.*r2.*y(3).*(1-p2) - d3.*y(4);
dydt = [tmp1; tmp2; tmp3; tmp4; p0];
end