Why won't Matlab plot or finish running my code due to the value of 1 variable?

I have a value known as h0 that I set as 10^-2 but whenever I run this with this value thge code never runs it just doesn'ty produce any values nor graphs. I am inputting this value as it represents mutants in this case y values that have a strong control over the population. Why wont this run?

this is what the code is based off of with tmp 1 being x0 and so on.

Code:

HematopoieticSystemModel

function [t,dydt1] = HematopoieticSystemModel

time = 0:3000;
%------%

x_0 = 1;
y_0 = 0;
x_1 = 0;
y_1 = 0;
x_2 = 0;
y_2 = 0;
% p0_0=100;
%----------%

p = [x_0 y_0 x_1 y_1 x_2 y_2];
%--------------------------------------------------%
[t,dydt1] = ode45(@myODE,time,p);


%---Matrices-----%
% [m1, m2, m3, m4, m5, m6] = equillibriatest(dydt1(:,1), dydt1(:,2), dydt1(:,3), dydt1(:,4), dydt1(:,5), dydt1(:,6));
% save('m1.mat','m1')
% save('m2.mat','m2')
% save('m3.mat','m3')
% save('m4.mat','m4')
% save('m5.mat','m5')
% save('m6.mat','m6')


%% -----Individual Plots------ %
% figure (1)
% hold on
% subplot (1,3,1);
% plot(t,dydt1(:, 1), 'r');
% title ('HSC');
% xlabel('time')
% ylabel('cells')
% hold off
% 
% hold on
% subplot (1,3,2);
% plot (t, dydt1(:,3),'g');
% title ('ST-HSCs');
% xlabel('time')
% ylabel('cells')
% hold off
% 
% hold on
% subplot (1, 3, 3);
% plot (t, dydt1(:,5),'b');
% title ('MPPs');
% xlabel('time')
% ylabel('cells')
% hold off

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

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

figure (3)
hold on
plot(t,dydt1(:, 2), 'g');
plot (t, dydt1(:,4),'b');
plot (t, dydt1(:,6),'m');
legend ('HSC' , 'ST-HSC', 'MPP');
xlabel('time (days)');
ylabel('mutant cells');
title ('Hematopoietic System: Mutant Growth');
set(gca, 'YScale', 'log');
hold off

figure(4)
hold on
plot (t, dydt1(:,1),'g');
plot (t, dydt1(:,2), 'k');
title ('HSCs');
xlabel('time');
ylabel('cells');
legend ('HSC', 'mutant');
set(gca, 'YScale', 'log');
hold off

figure (5)
hold on
plot (t, dydt1(:,3),'b')
plot (t, dydt1(:,4),'k')
title ('ST-HSCs');
xlabel('time');
ylabel('cells');
legend ('ST-HSC', 'mutant');
set(gca, 'YScale', 'log');
hold off

figure (6)
hold on
plot (t, dydt1(:,5), 'm')
plot (t, dydt1(:,6),'k')
title ('MPPs');
xlabel('time');
ylabel('cells');
legend ('MPPs', 'mutants');
set(gca, 'YScale', 'log');
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);

global P0 P1 P2 h0 r0 r1 r2 d3 u1 s x0 u0 u2 h2 h1 
d3=0.04;
s=0.03;
x0=17000;
r0=0.01818;
r1=0.8712;
r2=8.0217;
P0=0.5;
P1=0.4783;
P2=0.4987;
u0=10^-5;
u1=10^-5;
u2=10^-5;
%%strong control
h0= 10^-2;
h1= 10^-2;
h2= 10^-2;

% y(1) = x_0;
% y(2) = y_0;
% y(3) = x_1;
% y(4) = y_1;
% y(5) = x_2;
% y(6) = y_2;

c0=0.5+((h0*x0)/2);
c1=P1+((r0*x0*h1*P1)/(r1*(1-2*P1)));
c2=P2+((2*r0*x0*h2*P2*(1-P1))/(r2*(1-2*P1)*(1-2*P2)));

p0=c0./(1+h0.*(y(1)+y(2)));
p1=c1./(1+h1.*(y(3)+y(4)));
p2=c2./(1+h2.*(y(5)+y(6)));

p_0=P0.*(1+s);
p_1=P1.*(1+s);
p_2=P2.*(1+s);


tmp1 = r0.*y(1).*p0.*(1-u0)-r0.*y(1).*(1-p0);
tmp2 = r0.*y(1).*p0.*u0-r0.*y(2).*(2.*p_0-1); 
%%
tmp3 = r0.*y(1).*(1-p0).*(2-u0)+r1.*y(3).*p1.*(1-u1)-r1.*y(3).*(1-p1);
tmp4 = r0.*y(1).*(1-p0).*u0+2.*r0.*y(2).*(1-p_0)+r1.*y(3).*p1.*u1+r1.*y(4).*(2*p_1-1);
%%
tmp5 =r1.*y(3).*(1-p1).*(2-u1)+r2.*y(5).*p2.*(1-u2)-r2.*y(5).*(1-p2);
tmp6 =r1.*y(3).*(1-p1).*u1+2.*r1.*y(4).*(1-p_1)+r2.*y(5).*p2.*u2+r2.*y(6).*(2*p_2-1);

dydt = [tmp1; tmp2; tmp3; tmp4; tmp5; tmp6];

end

function [matrix1, matrix2, matrix3, matrix4, matrix5, matrix6] = equillibriatest (in1, in2, in3, in4, in5, in6)
matrix1 = in1;
matrix2 = in2;
matrix3 = in3;
matrix4 = in4;
matrix5 = in5;
matrix6 = in6;
end

Answer 1

[Walter Roberson]

You have to reduce the upper bound on the time by quite a bit. With you current code, with ode45, after running for about 45 minutes, it would not have proceeded much past 42 simulated seconds.

If you switch to ode15s, then it quits the integration after a small fraction of a second.

If you switch to ode23s, then after roughly 0.097 seconds it really starts to slow down; 0.098 as an upper bound is tolerable but beyond that it gets very slow.

With ode45 an upper bound of 20 seconds is tolerable time.

With both ode23s and ode45, dydt1(:,3) and dydt1(:,4) go negative early on, and so are mostly not plotted because you set log axes.

You also have a little bit of a settings problem. You hold on before drawing anything at all in the plots, and that locks in a ylim of [0 1]. You should use code closer to

figure (2)

plot(t,dydt1(:, 1), 'g');
hold on
plot (t, dydt1(:,3),'b');
plot (t, dydt1(:,5),'m');
legend ('HSC' , 'ST-HSC', 'MPP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off
ylim auto

[Walter Roberson]

I ran with ode45 up to time t = 20. It took a while. I used [0, 20] rather than 0:20 for the time, so it emitted an output (four outputs actually) each time it took a successful step. It emitted over 12 million output rows up to t = 20.

dydt1(:, 3) and dydt1(:, 4) both go negative pretty much immediately, and they are both quite unstable -- nearly all of the hard work is in correctly tracing :3 and :4 . The :3 is pretty much pure noise; the :4 is noise overlayed on a distinct arc.

I would suggest that you re-check that your implementation matches your equations.