a=2; b=0.5; l=10;
modell={x'[t]==a -(b+1) x[t] + x[t]^2 y[t],
y'[t]==b x[t] - x[t]^2 y[t],
x[0]==u0, y[0]==v0};
p:=Module[{},loesung=NDSolve[modell,{x,y},{t,0,l}]; 
abb=ParametricPlot[Evaluate[{x[t],y[t]}/.loesung],{t,0,l},
DisplayFunction->Identity,
PlotPoints->200]
];
max=5;
bild1a=Table[{u0=max i/max;v0=max;p},{i,0,max,1}];
bild1b=Table[{u0=max;v0=max j/max;p},{j,0,max,1}];
bild2a=Table[{u0=max i/max;v0=0.1;p},{i,0,max,1}];
bild2b=Table[{u0=0.1;v0=max j/max;p},{j,0,max,1}];
abb0=Show[Evaluate[bild1a,bild1b,bild2a,bild2b],
DisplayFunction->$DisplayFunction,
PlotRange->All,AxesLabel->{"x","y"},
ImageSize->{400,300}];