How to use 5 function coupled each other using ODE45? it is possible?

I have a coupled differential equation. I'm confused why my M3, O, and P values are 0. Even though the initial conditions M2, M3, O and P are 0 but only M2 has a value. Is there something wrong with the function?

function CM1 = mymode (t,M1,M2,M3,O,P)
M1= 10; 
M2 = 0;
M3 = 0; 
O  = 0; 
P=0;

delta=50;
gamma=75;
K1= 10^-4;
K2=5*10^-4;
K3=10^-3;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;

CM1= zeros(5,1);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((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)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);


end


[t,M1,M2,M3,O,P] = ode45(@mymode, [0,100],[0,0.01])
plot (t,M1,M2,M3,O,P)

[Torsten]

The list of inputs to "myode" is a scalar (t) and the vector of unknowns (y). Thus you have to pass your unknowns as a vector, not as single separated scalar values. @Alan Stevens showed you how to do this.

Further, by setting

function CM1 = mymode (t,M1,M2,M3,O,P)
M1= 10; 
M2 = 0;
M3 = 0; 
O  = 0; 
P=0;
...    

you overwrite the solution variables and calculate the derivatives of your equations for M1 = 10, M2 = 0, M3 = 0, O = 0 and P = 0 for all times t. I doubt this is what you want.

[cindyawati cindyawati]

The results that Alan did were what I wanted, but my lecturer tried to use order 4 runge kutta. I'm still stuck if use order 4 runge kutta, Is it possible that ODE45 can be applied in combination with order 4 runge kutta? @Torsten

[Torsten]

This is a Runge-Kutta-4 code for your problem. Try to understand how "runge_kutta_RK4" works on your system of equations to do better next time.

tstart = 0.0;
tend = 100.0;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0];
f = @myode;
Y = runge_kutta_RK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
O  = Y(:,4);
P  = Y(:,5);
figure
subplot(3,1,1)
plot(T,M1),grid, title('M1')
subplot(3,1,2)
plot(T,M2),grid, title('M2')
subplot(3,1,3)
plot(T,M3),grid, title('M3')

figure
subplot(2,1,1)
plot(T,O),grid, title('O')
subplot(2,1,2)
plot(T,P),grid, title('P')

function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-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(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end
function CM1 = myode (~,MM)
  M1 = MM(1);
  M2 = MM(2);
  M3 = MM(3);
  O  = MM(4);
  P  = MM(5);
  delta=50;
  gamma=75;
  K1= 10^-4;
  K2=5*10^-4;
  K3=10^-3;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  CM1= zeros(1,5);
  CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((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)-(Ko*M1*O)-(mu_o*O);
  CM1(5) = (Ko*M1*O)-(mu_p*P);
end

[cindyawati cindyawati]

thank you @Torsten

Answer 1

[Alan Stevens]

Better like this:

MM0 = [10, 0, 0, 0, 0];
tspan = [0 100];

[t, MM] = ode15s(@mymode, tspan,MM0);
M1 = MM(:,1); 
M2 = MM(:,2);
M3 = MM(:,3); 
O  = MM(:,4); 
P  = MM(:,5);
figure
subplot(3,1,1)
plot(t,M1),grid, title('M1')
subplot(3,1,2)
plot(t,M2),grid, title('M2')
subplot(3,1,3)
plot(t,M3),grid, title('M3')

figure
subplot(2,1,1)
plot(t,O),grid, title('O')
subplot(2,1,2)
plot(t,P),grid, title('P')

function CM1 = mymode (~,MM)
M1 = MM(1); 
M2 = MM(2);
M3 = MM(3); 
O  = MM(4); 
P  = MM(5);

delta=50;
gamma=75;
K1= 10^-4;
K2=5*10^-4;
K3=10^-3;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;

CM1= zeros(5,1);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((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)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);


end

[cindyawati cindyawati]

Thank you for your response @Alan Stevens This code use ode15? Can't use ode 45?

[Alan Stevens]

Yes, you can also use ode45 - it just uses more points over the time span.

[cindyawati cindyawati]

Oke thank you