This is the basic Lotka-Volterra model similar to the code that we assembled in InsightMaker.

The differential equations looks like

x’(t)=alpha x (1- x/(alpha x/(beta y)))

where (alpha x/(beta y)) is the “effective” carrying capacity of the prey.

y’(t)=delta y (x - gamma/delta)

is the predator (lynx) equation. 

Author: Andy Long, Northern Kentucky University (2020)

Clear["Global`*"]

Manipulate[
lv=
NDSolve[
{(* x is prey: *)
x'[t]==alpha  x[t](1-x[t]/((alpha*x[t])/(beta*y[t]))) ,
(* y is predator: *)
y'[t]==delta  y[t](x[t]- gamma/delta),
x[0]==gamma/delta*popinitial,
y[0]==alpha/beta*popinitial
},
{x[t], y[t]},
{t,0,tend}];
xeq=gamma/delta;
yeq=alpha/beta;
popvals=Table[Evaluate[{x[t]/xeq,y[t]/yeq}/.lv],{t,0,tend,.1}];
p1=Plot[Evaluate[{x[t]/xeq}/.lv],{t,0,tend},
PlotRange->All,AxesOrigin->{0,0},PlotLabel->"Populations (Lynx scaled by q) relative to K",
AxesLabel->{Years,Pop}];
p2=Plot[Evaluate[{y[t]/yeq}/.lv],{t,0,tend},PlotStyle->{Orange,Dashed,Thick}];
Show[GraphicsGrid[{{ListPlot[popvals,PlotLabel->"Phase Portrait, relative to equilibria",
AxesLabel->{"Hare","Lynx"}]},{Show[p1,p2]}}]],
{{lv,{}},None},
{{tend,28},2,100,1},
{{popinitial,1.01},.1,2,.1},
{{alpha,1},0.1,20,.1},
{{beta,0.035},.001,1,.1},
{{gamma,1},.01,2,.01},
{{delta,0.01},.001,1,.01}
]