For β =0.75

In[]:=
r=1;​​b=1;​​c=0.01;​​a=0.36;​​μ=0.4;
In[]:=
f[x_,y_]:=r*x-b*x^2-c*x*y-0.75*x*y/(a+x);​​g[x_,y_]:=-μ*y+0.75*x*y/(a+x);
In[]:=
Q[a_,b_,h_,N_]:=(u[0]=a;v[0]=b;​​Do[{u[n+1]=u[n]+h*f[u[n]+h/2*f[u[n],v[n]],v[n]+h/2*g[u[n],v[n]]],v[n+1]=v[n]+h*g[u[n]+h/2*f[u[n],v[n]],v[n]+h/2*g[u[n],v[n]]]},{n,0,N-1}]);
In[]:=
Q[0.2,0.05,1,200];
In[]:=
X=Interpolation[Table[{n,u[n]},{n,0,200}]];​​Y=Interpolation[Table[{n,v[n]},{n,0,200}]];
In[]:=
Plot[{X[t],Y[t]},{t,0,200}]
Out[]=
50
100
150
200
0.0
0.2
0.4
0.6
0.8
In[]:=
Clear["Global`*"];

For β = 1

In[]:=
r=1;​​b=1;​​c=0.01;​​a=0.36;​​μ=0.4;​​
f[x_,y_]:=r*x-b*x^2-c*x*y-1*x*y/(a+x);​​g[x_,y_]:=-μ*y+1*x*y/(a+x);​​
In[]:=
Q[a_,b_,h_,N_]:=(u[0]=a;v[0]=b;​​Do[{u[n+1]=u[n]+h*f[u[n]+h/2*f[u[n],v[n]],v[n]+h/2*g[u[n],v[n]]],v[n+1]=v[n]+h*g[u[n]+h/2*f[u[n],v[n]],v[n]+h/2*g[u[n],v[n]]]},{n,0,N-1}]);​​
In[]:=
Q[0.2,0.05,1,200];​​
In[]:=
X=Interpolation[Table[{n,u[n]},{n,0,200}]];​​Y=Interpolation[Table[{n,v[n]},{n,0,200}]];​​
In[]:=
Plot[{X[t],Y[t]},{t,0,200}]
Out[]=
50
100
150
200
0.2
0.4
0.6
0.8