Numerical Analysis of a Phase Portrait

Definition of a vector field

f[x_,y_]:=1-(b+1)*x+a*x^2*y;
g[x_,y_]:=b*x-a*x^2*y;

Calculation of Fixed Points

Simplify[Solve[{f[x,y]==0,g[x,y]==0},{x,y}]]

Calculation of Jacobian

Simplify[D[f[x,y],x]]
Simplify[D[f[x,y],y]]
Simplify[D[g[x,y],x]]
Simplify[D[g[x,y],y]]

J[x_,y_]:=(({{ -1-b+2 a x y, a (x^2)}, {b-2 a x y, -a (x^2)}}))

Calculation of Eigenvalues and stability of fixed points

Eigenvalues[J[1,b/a]]

Since we know that the characteristic polynomial is given by p(λ)=(λ^2)-Tr(A)λ+Det(A) it follows that 
(Subscript[λ, 1,2])=Tr(A)±(Sqrt[((Tr(A))^2)-4Det(A)])=Tr(A)±Tr(A)(Sqrt[1-4Det(A)/((Tr(A))^2)])
Therefore, a fixed point is stable if Tr(A)<0 and Det(A)>0 and is unstable otherwise. 
However, the eigenvalues can be complex leading to oscillations. Specifically, if 4Det(A)>((Tr(A))^2)the eigenvalues are complex. There are two cases
1. If Tr(A)<0 and 4Det(A)>((Tr(A))^2) the fixed point is classified as a stable spiral
2. If Tr(A)>0 and 4Det(A)>((Tr(A))^2) the fixed point is classified as an unstable spiral

This means we can solve for this conditions to determine different cases for the phase portraits

Tr[J[1,b/a]]

Det[J[1,b/a]]

Therefore the fixed point is stable if a>b-1 and a>0.

Moreover, the eigenvalues will be complex if 4a>(1+a-b)^2.

Case 1: Tr(A)<0, Det(A)>0 and 4Det(A)>((Tr(A))^2).

a=2;
b=2;
xlowlim=0;
xuplim=2;
ylowlim=0;
yuplim=2;
plot1=StreamPlot[{f[x,y],g[x,y]},{x,xlowlim,xuplim},{y,ylowlim,yuplim},PlotRange->{{xlowlim,xuplim},{ylowlim,yuplim}},FrameLabel->{"x","y"},LabelStyle->Directive[Black, Large]];
plot2=ContourPlot[f[x,y]==0,{x,xlowlim,xuplim},{y,ylowlim,yuplim},ContourStyle->{Thick,Blue}];
plot3=ContourPlot[g[x,y]==0,{x,xlowlim,xuplim},{y,ylowlim,yuplim},ContourStyle->{Thick,Red}];
s=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==0,y[0]==1/2},{x,y},{t,20}];
plot4=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20},PlotStyle->{Thick,Black}];
s=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==1/2,y[0]==0},{x,y},{t,20}];
plot5=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20},PlotStyle->{Thick,Black}];
Show[plot1,plot2,plot3,plot4,plot5]
Plot[{Evaluate[x[t]/.s],Evaluate[y[t]/.s]},{t,0,20},PlotRange->All,AxesLabel->{"t","x(t),y(t)"},LabelStyle->Directive[Black,14]]
Clear[a,b]

This is a stable spiral.

Case 2: Tr(A)>0, Det(A)>0 and 4Det(A)>((Tr(A))^2).

a=1;
b=2.5;
xlowlim=0;
xuplim=4;
ylowlim=0;
yuplim=4;
plot1=StreamPlot[{f[x,y],g[x,y]},{x,xlowlim,xuplim},{y,ylowlim,yuplim},PlotRange->{{xlowlim,xuplim},{ylowlim,yuplim}},FrameLabel->{"x","y"},LabelStyle->Directive[Black, Large]];
plot2=ContourPlot[f[x,y]==0,{x,xlowlim,xuplim},{y,ylowlim,yuplim},ContourStyle->{Thick,Blue}];
plot3=ContourPlot[g[x,y]==0,{x,xlowlim,xuplim},{y,ylowlim,yuplim},ContourStyle->{Thick,Red}];
s=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==0,y[0]==1/2},{x,y},{t,20}];
plot4=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20},PlotStyle->{Thick,Black}];
s=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==1/2,y[0]==0},{x,y},{t,20}];
plot5=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20},PlotStyle->{Thick,Black}];
Show[plot1,plot2,plot3,plot4,plot5]
Plot[{Evaluate[x[t]/.s],Evaluate[y[t]/.s]},{t,0,20},PlotRange->All,AxesLabel->{"t","x(t),y(t)"},LabelStyle->Directive[Black,14]]
Clear[a,b]

This is an unstable spiral surrounded by a limit cycle.