find equilirium point with all parameter known

Hello

I have a mathematical model for degree 7 non-linear system. I have calculated the equilibrium points for a disease-free system (x2=0 and x3=0).

However, it is not possible to find endemic equilibrium points manually, so it is necessary to use numerical methods. but I don't know how to use matlab's feature to get the endemic equilibrium point (x2≠0 , x3≠0) .

can you help me change this script, here i already have the parameter value.

Pg=0.0068;
Ps=0.012;
Pr=0.006;%0.012 0.002 0.006
pe=0.002;
beta1=0.00702;
beta2=0.00702;
beta3=0.00702;
d4=1.7143;
u=0.05;
a1=1.10;
a2=4.6205;
a4=4.6666;
a5=1.10;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
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.01813;
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 = 1000;

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,u,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma, alpha, Pr);

function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,u,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));

v=0.0;
Phi=3.3*10^-3;

    H=0;
if x(6)>0 || x(7)>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

Answer 1

[Torsten]

I don't know if this solves your problem, but of course you have to reverse the order of the commands here:

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));

You cannot use the if-statement using hdot(4) before it is computed.