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
