Final figures

Preparations: Parabolas

fig[A0_,A1_,A2_,A3_,from_,to_,number_,letter_]:=Module[{param={A0,A1,A2,A3}},H[α_]:=A3 (α^3)+A2 (α^2)+A1 α+A0;
Plot[H[α],{α,from,to},AspectRatio->1,AxesOrigin->{0,0},PlotLabel->Row[{number,letter}], Ticks->None,Epilog->{PointSize[Large],Red,Point[{-(A2/(3A3)),H[-(A2/(3A3))]}]},ImageSize->200]]

Fig. 1

GraphicsArray[{a1=fig[1+(1/24),11/6,1,1/6,-3.2,0.4,1,a],a2=fig[1,1,1,1/6,-6,2.5,1,b],a3=fig[2,2/3,1/4,1/12,-6,2.5,1,c]}]

}

Fig. 2

GraphicsArray[{a4=fig[1,5/12,-1/2,1/12,-1.5,5,2,a],a5=fig[1,9/16,-3/8,1/16,-1.8,5,2,b],a6=fig[1,4,-3/2,1/4,-1,4,2,c]}]

Fig. 3

GraphicsArray[{a7=fig[10/15,-16/15,3/15,1/15,-7,4,3,a],a8=fig[11/10,-9/10,1/10,1/10,-4.5,3,3,b]}]

Fig. 4

GraphicsArray[{a9=fig[12/10,-4/10,-3/10,1/10,-2.5,4,4,a],a10=fig[22/10,-4/10,-3/10,1/10,-2.5,4,4,b]}]

Preparations: Typical null clines of prey in case of two dimensions

sol=Solve[(r x (1-(x/K))-y m(x/(a y+x)))==0,y][[1,1]];

pl[R_,from_,opts___]:=Plot[y/.sol/.{r->R,K->0.1,m->16,a->4},{x,from,0.10},AxesLabel->{x,(Subscript[y, 1])},AspectRatio->1,opts]

p1=pl[4,0,Epilog->Text[Style["m/a=r",Large],{0.03,0.012}]];
p2=pl[7,0.0,PlotRange->{0,0.15},Epilog->Text[Style["m/a<r",Large],{0.03,0.12}]];
p3=pl[3.5,0,Epilog->Text[Style["m/a>r",Large],{0.03,0.0045}]];

Fig. 5: Typical null clines

Show[GraphicsArray[{p1,p2,p3}],ImageSize->600,PlotLabel->"Typical null clines"]

Preparations: Typical null clines of prey in case of three dimensions

ClearAll[param,sh,stac,r,K,m1,a1,m2,a2];
param[R_]:={r->R,K->0.1`,m1->16,a1->4,m2->18,a2->2};
stac=Solve[{r (1-(x/K))-y1 m1(1/(a1 y1+x))-y2 m2(1/(a2 y2+x))==0,m1(x/(a1 y1+x))-d1 ==0,m2(x/(a2 y2+x))-d2==0},{x,y1,y2}];sh[R_,xig_,y1ig_,y2ig_,plp_:2]:=Show[ContourPlot3D[((r x (1-(x/K))-y1 m1(x/(a1 y1+x))-y2 m2(x/(a2 y2+x)))/.param[R])==0,{x,0,xig},{y1,0,y1ig},{y2,0,y2ig},AxesLabel->{x,(Subscript[y, 1]),(Subscript[y, 2])},PlotPoints->plp],Graphics3D[{RGBColor[0,0,0],Sphere[({x,y1,y2}/.stac/.param[R]/.{d1->8,d2->12})[[1]],0.002]}],PlotLabel:>ToString[param[R]]<>"\n",LabelStyle->Directive[Bold]]

Fig. 6

sh[3,0.1,0.01,0.005]

Fig. 7

sh[7,0.1,0.035,0.025]

Fig. 8

sh[10,0.1,0.035,0.08]

Summarizing (Be patient!)

Manipulate[sh[R,0.1,0.035,0.05],{{R,7},3,14,1}]

Preparations: Time evolution of species

ClearAll[rhs,par,stac,r,α];
rhsreduced={r (1-(x/K))-y1 m1(1/(a1 y1+x))-y2 m2(1/(a2 y2+x)),m1(q/(a1 y1+q))-d1,m2(q/(a2 y2+q))-d2,α(x-q)};
par={K->1/10,m1->16,a1->4,m2->18,a2->2,d1->8,d2->12};
stac=Solve[Thread[(rhsreduced/.par)==0],{x,y1,y2,q}][[1]]

figs[r_,α_,x0_,y10_,y20_,q0_,from_,to_]:=Module[{nds},
Needs["PlotLegends`"];
nds=NDSolve[{x'[t]==x[t](r(1-((x[t])/K))-y1[t] m1(1/(a1 y1[t]+x[t]))-y2[t] m2(1/(a2 y2[t]+x[t]))),y1'[t]==y1[t](m1((q[t])/(a1 y1[t]+q[t]))-d1),y2'[t]==y2[t](m2((q[t])/(a2 y2[t]+q[t]))-d2),q'[t]==α(x[t]-q[t]),x[0]==x0,y1[0]==y10,y2[0]==y20,q[0]==q0}/.par,{x,y1,y2,q},{t,0,100}];px=Plot[Tooltip[x[t]/.nds,"x[t]"],{t,from,to},PlotStyle->{Directive[Red]},PlotRange->All,AxesOrigin->{0,0},ImageSize->300];
py1=Plot[Tooltip[y1[t]/.nds,"(Subscript[y, 1])[t]"],{t,from,to},PlotStyle->{Directive[Yellow,Thick]},PlotRange->All,AxesOrigin->{0,0},ImageSize->300];
py2=Plot[Tooltip[y2[t]/.nds,"(Subscript[y, 2])[t]"],{t,from,to},PlotStyle->{Directive[Blue,Dashed]},PlotRange->All,AxesOrigin->{0,0},ImageSize->300];
pq=Plot[Tooltip[q[t]/.nds,"q[t]"],{t,from,to},PlotStyle->{Directive[Green,Thick]},PlotRange->All,AxesOrigin->{0,0},ImageSize->300];
GraphicsArray[{Show[px,py1,py2,pq],Show[ParametricPlot3D[{x[t],y1[t],y2[t]}/.nds,{t,from,to},BoxRatios->{1,1,1},PlotRange->All,AxesLabel->{"\nx","\!\(\*SubscriptBox[\(y\), \(\(1\)\(\\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \)\)]\)","(Subscript[y, 2])                  "},TicksStyle->Directive[Orange,10],LabelStyle->Directive[Blue,Medium]],
Graphics3D[{Red,PointSize[Large],Point[{((-5+r)/(10 r)),((-5+r)/(40 r)),((-5+r)/(40 r))}]}],ImageSize->300]}]]

Preparation: H(α) function

rhs=rhsreduced {x,y1,y2,1};

rhs//MatrixForm

stac

(jacobian=D[rhs,{{x,y1,y2,q}}]//Simplify)//MatrixForm

(matrix=jacobian/.stac/.par//Simplify)//MatrixForm

cp[λ_]:=CharacteristicPolynomial[matrix,λ]

cp[λ]

{A0,A1,A2,A3,A4}=CoefficientList[cp[λ],λ]

ClearAll[H];
H[α_]:=A3(A1 A2-A0 A3)-(A1^2)

Collect[H[α],α]

H[α]/.r->13

ClearAll[inflexion];
inflexion[r_]:=Solve[D[H[α],{α,2}]==0,α][[1,1,2]]

Point[{inflexion[13],0}]

{inflexion[13]/.{r->13},H[α]/.{α->inflexion[13]}/.r->13}

Fig. 9

figs[13,1,0.1,0.1,0.1,0.1,0,10]

Fig. 10

Plot[Evaluate[H[α]/.{r->13}],{α,-6,1},Epilog->{PointSize[Medium],Point[{inflexion[13]/.{r->13},H[α]/.{α->inflexion[13]}/.r->13}]}]

Fig. 11

figs[7,1,(3/100),(7/1000),(7/1000),(2/100),0,10]

Fig. 12

Plot[Evaluate[H[α]/.{r->7}],{α,-5,4},Epilog->{PointSize[Medium],Point[{inflexion[7]/.{r->7},H[α]/.{α->inflexion[7]}/.r->7}]}]

Fig. 13

figs[7,1,(3/100),(7/1000),(7/1000),(2/100),0,70]

Fig. 14

figs[8,0.2,0.03,0.009,0.009,0.03,0,23]

Fig. 15

Manipulate[Module[{nds,x,y1,y2,q},nds=NDSolve[{x'[t]==x[t](r(1-((x[t])/0.1))- ((16y1[t])/(4 y1[t]+x[t]))- ((18y2[t])/(2 y2[t]+x[t]))),y1'[t]==y1[t](((16q[t])/(4 y1[t]+q[t]))-8),y2'[t]==y2[t](((18q[t])/(2 y2[t]+q[t]))-12),q'[t]==α(x[t]-q[t]),x[0]==0.03,y1[0]==0.009,y2[0]==0.009,q[0]==0.03},{x,y1,y2,q},{t,0,50},MaxSteps->10000][[1]];
ParametricPlot[Evaluate[{x[t],y1[t]}/.nds],{t,0,100},PlotStyle->{Red,Thickness[0.01]},AspectRatio->1,AxesLabel->{"x","(Subscript[y, 1])"}, PlotRange->{{0,0.1},{0.006,0.015}},Epilog->{PointSize[Large],Blue,Point[{0.03,0.009}]}]],{{α,0.25},0.1,0.3,0.05},{{r,8},7,9,0.1},LocalizeVariables->False]