Solving system with ode45

%I'm trying to solve the following system for an array with 10 elements (k=10) for biology project:

% Most of the code includes paramters and while loops, but I'm wondering
% Is it fine? How to create figures for each of the variables or together? 
% I need to get 10 lines for each variable. Thanks




etana=0.001;
etama=1;
etaag=0.2;
etagm=1;
etaap=0.1;
ron=0;
roa=1.5;
rog=0.2;
rom=0;
rop=0;
mun=0.5;
mua=0.2;
mug=0.1;
mum=0.01;
mup=0.3;
beta=20;
Bmax=10^6;
f=10;
gamma=0.08;
delta=0.09;
tspan = [0 100];
y0 = [100 0 0 0 0 1 0 ];
si=1;
u=0;

while u<101
    u=u+100;
    tspan=[u-100 u];
   
[t,y] = ode45(@(t,y) odefcn(t,y,etana,etama,etaag,etagm,etaap,ron,roa,rog,rom,rop,mun,mua,mug,mum,mup,beta,Bmax,f,si,gamma,delta), tspan, y0);
y0=y(end,:);
y0(6)=1;

plot(t,y(:,6),'-',t,y(:,7),'.-')
hold on
end
hold off 

function dydt = odefcn(t,y,etana,etama,etaag,etagm,etaap,ron,roa,rog,rom,rop,mun,mua,mug,mum,mup,beta,Bmax,f,si,gamma,delta)
  dydt = zeros(7,1);
  N=zeros(10,1);
  A=zeros(10,1);
   G=zeros(10,1);
  M=zeros(10,1);
  P=zeros(10,1);
  k=1;
  while k<=10
       si=1/(k^2);
  N(k)=y(1);
  A(k)=y(2);
  G(k)=y(3);
  M(k)=y(4);
  P(k)=y(5);
  R=y(6);
  RG=y(7);
  dydt(1) =beta+(ron-mun-R*si*etana)*N(k) ;
  dydt(2) =R*si*etana*N(k)+R*si*etama*M(k)+(R*roa-mua-R*si*(etaag+etaap))*A(k);
  dydt(3) =R*si*etaag*A(k)+(rog*si*RG-mug-etagm)*G(k);
  dydt(4) =etagm*G(k)+(rom-mum-R*si*etama)*M(k);
  dydt(5) =R*si*etaap*A(k)+(rop-mup)*P(k);
  dydt(6) =dydt(6)-R*(A(k)+f*P(k))*si/Bmax;
 dydt(7)=gamma*R-delta*RG;
 k=k+1;
  end
end

Answer 1

[Star Strider]

‘I need to get 10 lines for each variable.’

This implies that one (or more) parameters must be varied 10 times and then integrated for each new set of parameters. The way to do that is to define the variables that need to change, then use a for loop to change them, and then run the integration in each iteration of the loop, with different parameter values (or combinations of values,depending on what the assignment requires). It would be easiest to save the outputs (‘y’ here) in a cell array to make the indexing easier. Then separate out the individual results after all the integrations are complete.

This will be easier if ‘tspan’ is a vector of more than two elements, for example —

tspan = linspace(0, 100, 50);

This sets all the outputs to the same length and the same times, making plotting them against the same time vector easier.

I would use one subplot for each variable, setting the xlim and ylim properties the same for all the subplot plots.

.