@aureux A and B are given by params as found above (now with g=0) . Thus if the original parameters a,b,c, ... are given, then A and B are determined. But during an investigation of the properties of the system (with a,b,c,.. not given) you can give them any value consistent with the assumptions made in the change of variables: a,b,k,beta all positive. The assumption c>=0 has nothing to do with that, but the wolves are sure going to die out if they don't profit from eating elk.
That implies A>=0 and B>0.

@aureux That looks good to me.

The advantage of making this change of variables from x,y,t to u,v,tau is that examining a system with only 2 parameters is easier than examining one with (now) 5 parameters.
The point being that the behaviors of sys and SYS are the same. Thus if you find that SYS has 3 equilibrium points, so does sys, and you can find the ones from the others.
The stability properties of those equilibrium points are the same: If the equilibrium point in the open first quadrant for SYS is asymptotically stable and a spiral point, so is the corresponding point for sys, etc.
After all (I guess) the question is: is the grey wolf - elk 'situation' in Yellowstone Park going to end in both species dying out, one of the two dying out, or can they coexist (not peacefully of course)?

@aureux Often people introduce dimensionless variables mostly because this usually will reduce the number of parameter (constants).
There is no unique way of doing this. I indicate a method below that uses dchange from the PDEtools package.
I replaced alpha with a, but otherwise didn't change names of constant in sys below.

The idea is make the substitutions x(t)= X*u(tau),y(t)=Y*v(tau), t=T*tau  thereby introducing also a new time tau. The 3 constants X,Y,T should be chosen so as to reduce the number of parameters. After all you have a,b,c,k,g,beta, i.e. 6 constants.
You can reduce that number to 2.
Suppose all your 6 parameters are nonnegative by assumption and suppose further that a>0, b>0, k>0, beta+g > 0, then the changes actually made below make sense. By that I mean that you don't want X,Y, or T to be zero (or negative even).

sys:=D(x)(t)=a*x(t)-a*x(t)^2/k-b*x(t)*y(t)-g*x(t), D(y)(t) = -beta*y(t)+c*x(t)*y(t)-g*y(t);

PDEtools:-dchange({x(t)=X*u(tau),y(t)=Y*v(tau),t=T*tau},{sys},[tau,u,v]);

solve(%,{diff(u(tau),tau),diff(v(tau),tau)});
SYS1:=collect(%,[u,v],distributed,factor);
## Now you can choose the constants X,Y,T the way you wish
## I did as follows
eval(SYS1,X=k/(a*T)); #Here you use that a>0 and k>0
eval(%,T=1/(beta+g)); #Here you use that beta+g>0
eval(%,Y=(beta+g)/b); #Here you use that b>0 and beta+g>0
eval(%,k=(a/c)*A); #Actually just introducing a simplifying parameter A
SYS:=op(eval(%,a-g=B*(beta+g))); # Introducing a simplifying parameter B
#The replacements made:
eqs:={X=k/(a*T),T=1/(beta+g),Y=(beta+g)/b,k=(a/c)*A,a-g=B*(beta+g)};
## To be able to get back and forth between sys and SYS we need these:
params:=solve(eqs,{T,X,Y,A,B});

From our assumptions it follows that A>=0, but B is just real.

@Yiannis Galidakis What do you want to happen here?

plot3d(x^2+y^2,x=-1..1,y=-1..1); p1:=%:
plot3d(x^2+y^2,x=-2..2,y=-2..2); p2:=%:
plots:-display(p1,p2,insequence);

########

Something like this?

f:=(x,y)->x^2+y^2;
plot3d(f,-1..1,-1..1); p1:=%:
plot3d((x,y)->f(x*2,y*2),-1..1,-1..1); p2:=%:
plots:-display(p1,p2,insequence);

@THAPELO There are a few misspellings and one syntax error:

(1) phi1 is spelled as ph1 in params.

(2) Sys1, Sys2 etc. are spelled sys1,sy2 in dsolve.

(4) a multiplication sign is missing in Sys5 right after theta.

Besides that you obviously have to include the information abut the meaning of N(t). This can be done simply by including it in params:

params:={Lambda1=0.029,b=0.5,phi1=0.25,phi2=1.01,mu=0.02,rho=0.1,delta=0.1,
f=0.85,theta=0.2,nu=0.33,h0=0.4,h1=0.08,h2=0.02,h3=0.06,h4=0.4,
  N(t)=S(t)+I1(t)+I2(t)+I3(t)+I4(t)+A(t)};

After that it works for me.

@Carl Love I believe you are wrong, see the example below. I remember reading about the change from rk4 to rkf45 some years ago and then thought the same as you are expressing.
However, I soon discovered that it does matter:

restart;
p1:=DEtools[DEplot](diff(x(t),t)=sin(x(t))+sin(t),x(t),t=0..100,[x(0)=0],linecolor=blue,arrows=none):
plots:-display(p1,size=[1000,200]);
p2:=DEtools[DEplot](diff(x(t),t)=sin(x(t))+sin(t),x(t),t=0..100,[x(0)=0],linecolor=blue,arrows=none,stepsize=0.1):
plots:-display(p2,size=[1000,200]);

Since the help page for DEplot doesn't mention the effect of stepsize even in the default case (rkf45) I shall suggest (via an SCR) that the help page be updated.

To see the limit cycle mentioned by me under the discussion using DEplot you don't need all these initial points.

Using one suffices. I chose one not too far away from the 3rd equilibrium point:

sol:=solve({eq1=0,eq2=0,eq3=0,eq4=0,eq5=0},{x,y,z,v,w});
pts:=map(subs,[sol],[x,y,z,v,w]);
eps:=.1:
ini:=[x,y,z,v,w](0)=~pts[3]+[eps$5];
res:=dsolve({Eq1,Eq2,Eq3,Eq4,Eq5} union {ini[]},numeric,maxfun=0);
plots:-odeplot(res,[x(t),y(t)],500..511);
XYZ:=plots:-odeplot(res,[x(t),y(t),z(t)],500..511):
YZV:=plots:-odeplot(res,[y(t),z(t),v(t)],500..511):
ZVW:=plots:-odeplot(res,[z(t),v(t),w(t)],500..511):
plots:-display(Array([XYZ,YZV,ZVW]));

@J4James I don't know what your objective is, but often people are interested in the long term behavior of the solutions (example: the Lorenz strange attractor).

If that is also your interest then you should use a time range like t=t1..t2, where t1 is somewhat large and t2-t1 is large enough to give an idea of the long term behavior. In case of a limit cycle that time difference should be just large enough to represent the cycle. In case of a more complicated long term behavior it is not so clear.

But start by trying  t = 250..260 or similar. While experimenting you might as well use stepsize=0.1.
Replace the default thickness with the one suggested by Carl.

There is graphical evidence of a limit cycle. To add more excitement I included z in the plot and used t=500..510.
Notice the need for increasing maxfun (actually I set it at infinity by maxfun=0):

DEtools[DEplot3d]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=500..510,initialset,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y,z],x=5.2..7.2,y=0.2..1.6,z=0..0.15,obsrange=false,thickness=0,maxfun=0);

@J4James I now remember that there is an option ,obsrange=false.

Use that and keep your original range.

@Carl Love I think you are right. I tried splitting the initial set in two according to the following procedure Q.
Notice that the range for x has been widened to begin with 5.

Q:=proc(ini) local x1,y1,a,b,c,d;
    a:=5; b:=7.2; c:=0.2; d:=1.6;
    x1,y1:=op(subs(ini,[x(0),y(0)])); (x1-a)*(x1-b)<=0 and (y1-c)*(y1-d)<=0
end proc;
Q(initialset[166]); #Test

iniset1,iniset2:=selectremove(Q,initialset):
nops~([iniset1,iniset2]);

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,x=5..7.2,y=0.2..1.6,iniset1 ,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,x=5..7.2,y=0.2..1.6,iniset2 ,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

@J4James I tried the following, where I only took the last 11 of the initial conditions and only used time 0..30 (my patience is limited too). I removed restrictions on x and y and also the arrows option as it doesn't make much sense here (and is ignored anyway).

It worked fine.

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,initialset[170..180],stepsize=0.01,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

@Vasudeva I see your point.

So if we let alpha = diff(beta(t),t) we get

ode:=diff(alpha(beta),beta)=1/alpha(beta)*sin(beta)*(r*alpha(beta)^2+g)/(2*r*cos(beta)-J);

#Using dsolve with the 'implicit' option, we get

dsolve(ode,implicit);

#and without:
sol:=dsolve(ode);

#Thus we have for the first one

ode2:=subs({alpha(beta)=diff(beta(t),t),beta=beta(t)},sol[1]);
dsolve(ode2);
res:=value(%);
odetest(res,eq);
#The second:
ode2b:=subs({alpha(beta)=diff(beta(t),t),beta=beta(t)},sol[2]);
dsolve(ode2b);
resb:=value(%);
odetest(resb,eq);

Without trying to understand what is going on I notice that you are using equality signs where assignments should be used.

Example:  You have x=x+1; That surely should be x:=x+1;
There are many such errors.

But there are other errors.

How did you do this by hand? By separation of variables? How?

@THAPELO You need to give the parameters to dsolve in one way or the other:

1. By assigning to the parameters before using dsolve, as in b:=1;  etc.
Be aware though that 'gamma' is Euler's constant in Maple, so you won't be allowed to assign to gamma. You could just call it something else, like gamma1.

2. A better way (in my opinion) is to define a set of equations (not assignments):

params:={b=1, d=0.2, k=1, alpha=0.5, gamma=0.25, mu=0.3};

and then use   eval({sys1, sys2, sys3},params)  as your concrete system. This is what I do below.
Using gamma is OK here as you are not assigning to it.

3. If you experiment a lot with varying parameters you may want to use the option parameters in dsolve

dsolve( {...ivp...}, numeric, parameters=[b, d, k, alpha, gamma, mu]);

This needs some explaining, which I shall not do now, but can give you an example if you wish.

##
If some parameters are fixed during the whole session you could assign to those as in point 1 and then use method 2 or 3 for the rest.

##############
restart;
sys1:= diff(S1(t),t) = b - d*S1(t) - k*S1(t)*S2(t)/(1+ alpha*S1(t)) - gamma*S3(t);
sys2:= diff(S2(t),t) = k*S1(t)*S2(t)/(1+ alpha*S1(t)) - (d + mu)*S2(t);
sys3:= diff(S3(t),t) = mu*S2(t) - (d + gamma)*S3(t);
##Method 2:
params:={b=1, d=0.2, k=1, alpha=0.5, gamma=0.25, mu=0.3};

Sol1:=dsolve(eval({sys1, sys2, sys3},params) union {S1(0)=1.35, S2(0)=0.9, S3(0)=0.45}, numeric);

plots:-odeplot(Sol1,[[t,S1(t)],[t,S2(t)],[t,S3(t)]],0..25);
#Values at t=25:
Sol1(25);
##Equilibria:
solve(rhs~({sys1, sys2, sys3}),{S1(t),S2(t),S3(t)});
eval([%],params);
map(subs,%,[S1(t),S2(t),S3(t)]);