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

a=-.01;
b=1;
xmax=4;
ymax=3;
x0=1.5;
y0=1.5;
Tmax=100;
plot1=ContourPlot[f[x,y]==0,{x,0,xmax},{y,0,ymax},FrameLabel->{"x","y"},FrameStyle->Directive[Black,Large],PlotRange->{{0,xmax},{0,ymax}},ContourStyle->{Thick,Red},PlotPoints->100];
plot2=ContourPlot[g[x,y]==0,{x,0,xmax},{y,0,ymax},ContourStyle->{Thick,Blue},PlotPoints->100];
plot3=StreamPlot[{f[x,y],g[x,y]},{x,0,xmax},{y,0,ymax}];
s=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==x0,y[0]==y0},{x,y},{t,Tmax}];
plot4=ParametricPlot[Evaluate[{x[t],y[t]}/. s],{t,0,Tmax},PlotStyle->{Thick,Black}];
Show[plot1,plot2,plot3,plot4]
Clear[a,b,xmax,ymax,x0,y0,Tmax];