trouble in plotting file differential equation contain function

I want to plot my script, glioma_if but it doesn't work, can you help me to plot?

close all; clear all; clc;

Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
Phi=2.0;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;

%t = linspace(0,0.1,100)';
%tspan = [0 1000];

X0 = [1.0 0.01 0.0 1.0 0.0];
t0 = 0;
tf = 800;

options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);

Index exceeds the number of array elements. Index must not exceed 5.

Error in solution>glioma0 (line 95)
if (x(6)>0)

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);


figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);

%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);

%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);

%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);

%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);

%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');

function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.0;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.001;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.01;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

[Torsten]

If you have 7 differential equations, you have to supply 7 initial conditions, not 5.

Answer 1

[Latifah]

I've added the initial condition, but it still error , can you fix it?

close all; clear all; clc;

Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
Phi=2.0;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;

%t = linspace(0,0.1,100)';
%tspan = [0 1000];

X0 = [1.0 0.01 0.0 1.0 0.0 0.01 0.01];
t0 = 0;
tf = 800;

options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);


figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);

%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);

%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);

%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);

%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);

%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');

function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.0;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.001;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.01;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

[Torsten]

@Latifah comment moved here:

Error using alpha

Too many output arguments.

Error in glioma_if>glioma0 (line 115)

hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in glioma_if (line 45)

[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5,

c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);

[Torsten]

Do you see a variable with name "alpha" somewhere in your code ? I cannot find it.

And this cannot be correct:

if (hdot(4)<0)
    hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
    hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

You use hdot(4) in an if-statement to compute hdot(3) although hdot(4) is not yet computed.

[Latifah]

finally got it to work thankyou, but I think it's still wrong. What if I just swap positions? can it be true? can you help me so that it can be calculated correctly?

close all; clear all; clc;

Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;
alpha=2.0;

%t = linspace(0,0.1,100)';
%tspan = [0 1000];

X0 = [1.0 0.01 0.0 1.0 0.0 0.01 0.01];
t0 = 0;
tf = 800;

options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma, alpha, Pr);
[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, psi, delta, gamma, alpha, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);


figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);

%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);

%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);

%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);

%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);

%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');

function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.0;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.001;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));

u=0.01;
Phi=200.0;

    H=0;
if (x(6)>0)
    H = 1;
%else%if (x(5)<=0);
%   H = 0;
end

hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));

if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);    
    
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end

hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));

if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);    
    
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end

hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);

hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1]; 
% To make xdot a column
% End of FUN1.M
end

[Torsten]

What if I just swap positions? can it be true?

Yes, of course. hdot(4) must be calculated before you use it in the if-statement to derive hdot(3).

But I would really be surprised if the integrator wouldn't come into trouble with all your if-statements in the function routines.

[Latifah]

okay thank you so much