modell = {s'[t] == - a s[t] i[t], i'[t]==a s[t] i[t] - b i[t], r'[t] == b i[t], s[0]==s0, i[0]==i0, r[0]==0}; a=0.0002; b=0.05; s0=250; i0=250; loesung = NDSolve[modell,{s,i,r},{t,0,100}]; abb=ParametricPlot[Evaluate[ {t,i[t]} /. loesung],{t,0,100}, AxesLabel->{"t","i"},ImageSize->{400,300}];