I have a mathematical model and its difficult to find endemic equilibrium point . can MATLAB do this

function hdot = glioma0(t,x,Pg, Ps, beta1, beta2, i1, i2, i3, a1, a2, a3, beta3, alpha, zeta, Pr);

% The time-dependent term is A * sin(w0 * t - theta)

%H = 1*sin(-1*t);

%G = x(1)

hdot=zeros(5,1);

%H= heaviside(x(5));

u=0.0;

Phi=200.0;

H=0;

if (x(5)>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))-(i1*x(1)*x(5))/(a1+x(1));

hdot(2)=Ps*x(2)*(1-x(2)-x(3))-beta2*x(1)*x(2)-u*x(2)*H-(i2*x(2)*x(5))/(a2+x(2));

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

if (hdot(1)<0)

hdot(4)=alpha*hdot(1)*x(4)-(i3*x(4)*x(5))/(a3+x(4));

elseif (hdot(1)>0)

hdot(4)=-(i3*x(4)*x(5))/(a3+x(4));

end

hdot(5)=Phi - zeta*x(5);

%hdot = hdot';

%hdot=[hdot(:);H,h1];

% To make xdot a column

% End of FUN1.M

end

this is the differential equation