dgl1=s'[t]==-β s[t] i[t]+f r[t]
dgl2=i'[t]==β s[t] i[t]-γ i[t] 
dgl3=r'[t]==γ i[t]-f r[t]

f=.1
β=.2
γ=.1

sol=NDSolve[{dgl1,dgl2,dgl3,s[0]==0.7,i[0]==0.3,r[0]==0},{s[t],i[t],r[t]},{t,0,1000}]

Manipulate[sol=NDSolve[
{
s'[t]==-β s[t] i[t] +f r[t],
i'[t]==β s[t] i[t] -γ i[t],
r'[t]==γ i[t] -f r[t],
s[0]==0.9,
i[0]==0.1,
r[0]==0
},
{s[t],i[t],r[t]},
{t,0,1000}];
Plot[
{
Evaluate[s[t]/.sol],
Evaluate[i[t]/.sol],
Evaluate[r[t]/.sol]
},{t,0,1000},PlotRange->All],
{f,0,.1},{β,0,.1},{γ,0,.1}]

p2=Plot[Evaluate[i[t]/.sol],{t,0,1000},PlotRange->All]

p3=Plot[Evaluate[r[t]/.sol],{t,0,1000},PlotRange->All]

Show[p1,p2,p3]

f1=0.0024
β1=0.0174
γ1=0.0036
sol1=NDSolve[
{
s'[t]==-β1 s[t] i[t] +f1 r[t],
i'[t]==β1 s[t] i[t] -γ1 i[t],
r'[t]==γ1 i[t] -f1 r[t],
s[0]==0.9,
i[0]==0.1,
r[0]==0
},{s[t],i[t],r[t]},
{t,0,1000}]

ParametricPlot[Evaluate[{s[t],i[t]}/.sol1],{t,0,1000}]

StreamPlot[{-β1 x y ,β1 x y -γ1 y},{x,0,1},{y,0,1}]

Manipulate[StreamPlot[{-β1 x y+f1 r ,β1 x y -γ1 y},{x,0,1},{y,0,1}],
{r,0,1}]