I have mathematical model with 7 equation . and its difficult to find endemic equilibrium point with analitical , i want to solve with numerical method

can you help me to find the endemic equilibrium point , with all the parameter known

close all; clear all; clc;

Pg=0.0068;
Ps=0.012;
Pr=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,y1]=ode45(@glioma1,[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 = glioma1(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.6;
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));

hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/(1+tau*x(4)))-beta3*x(1)*x(3)+u*x(2)*H-v*H*x(3);    

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

[Walter Roberson]

The mathematics of ode45 is not valid when you have discontinuities in the first or second deriviatives of the equations. You need to construct event functions to detect each condition under which you are changing your equations -- detect x(6) or x(7) changing between negative and positive, detect hdot(1) changing between negative and positive.

Question: what should your hdot(5) be if hdot(1) is exactly 0?