l=100; n1=200; n2=150; r1=0.00025; r2=0.00021; a1=0.05 ;a2=0.049;
modell={u'[s]==r1 v[s] (n1 - u[s]) - a1 u[s],
v'[s]==r2 u[s](n2 -v[s]) - a2 v[s] ,u[0]==u0, v[0]==v0};
p:=Module[{},loesung=NDSolve[modell,{u,v},{s,0,l}]; 
   abb=ParametricPlot[Evaluate[{u[s],v[s]}/.loesung],{s,0,l},
   DisplayFunction->Identity,
   PlotPoints->200]
];
max=150;
bild1a=Table[{u0=max i/max;v0=max;p},{i,1,max,25}];
bild1b=Table[{u0=max;v0=max j/max;p},{j,1,max,25}];
bild2a=Table[{u0=max i/max;v0=1;p},{i,1,max,25}];
bild2b=Table[{u0=1;v0=max j/max;p},{j,1,max,25}];
abb0=Show[Evaluate[bild1a,bild1b,bild2a,bild2b],
DisplayFunction->$DisplayFunction,
PlotRange->All,AxesLabel->{"I","I*"},ImageSize->{400,300}];