Never make changes in the RK4 function. Changes have to be made in the call to RK4 and in f.
tstart = 0;
tend = 100;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0 0 0 0];
f = @myode;
Y = fRK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
M4 = Y(:,4);
M5 = Y(:,5);
M6 = Y(:,6);
O = Y(:,7);
P = Y(:,8);
figure
%subplot(3,1,1)
%hold on
plot(T,M1,'r','Linewidth',2)
xlabel('Times (days)')
ylabel('M1 (gr/ml)')
figure
%subplot(3,1,2)
%hold on
plot(T,M2,'b','Linewidth',2)
xlabel('Times (days)')
ylabel('M2 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M3,'g','Linewidth',2)
xlabel('Times (days)')
ylabel('M3 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M4,'b--','Linewidth',2)
xlabel('Times (days)')
ylabel('M4 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M5,'r--','Linewidth',2)
xlabel('Times (days)')
ylabel('M5 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M6,'g--','Linewidth',2)
xlabel('times (days)')
ylabel('M6 (gr/ml)')
figure
%subplot(2,1,1)
%hold on
plot(T,O,'k','Linewidth',2)
xlabel('Times (days)')
ylabel('O (gr/ml)')
figure
%subplot(2,1,2)
%hold on
plot(T,P,'m','Linewidth',2)
xlabel('Times (days)')
ylabel('P (gr/ml)')
function Y = fRK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for j = 2:N
t = T(j-1);
y = Y(j-1,:);
h = T(j) - T(j-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(j,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
function CM1 = myode (~,MM)
M1 = MM(1);
M2 = MM(2);
M3 = MM(3);
M4 = MM(4);
M5 = MM(5);
M6 = MM(6);
O = MM(7);
P = MM(8);
delta=50;
gamma=75;
K1=10^-4;
K2=5*10^-4;
K3=10^-3;
K4=5*10^-3;
K5=10^-2;
K6=5*10^-2;
Ko=0.1;
n=6;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_4=10^-3;
mu_5=10^-3;
mu_6=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
sumterm = K2*M2 + K3*M3 + K4*M4 + K5*M5;
CM1= zeros(1,8);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*sumterm-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
CM1(4) = (K3*M1*M3)-(K4*M1*M4)-(mu_4*M4);
CM1(5) = (K4*M1*M4)-(K5*M1*M5)-(mu_5*M5);
CM1(6) = (K5*M1*M5)-(K6*M1*M6)-(mu_6*M6);
CM1(7) = (K6*M1*M6)-(Ko*M1*O)-(mu_o*O);
CM1(8) = (Ko*M1*O)-(mu_p*P);
end
