Clear[b, g, v, S, II, R]
  

b  =  
        0.00003;
g =  0.1587;
v =  
        0.034;
pi = 
        100.0;
u =  .0125;

ode1  =  S'[t] == 
        pi  -  u S[t] -  b S[t] II[t] +  
          g  R[t]
ode2  =  II'[t] == 
        b S[t] II[t] - v II[t]  -  u II[t]
ode3  =  R'[
          t] == v II[t] - g R[t]

Solve[{- b S II +  g  R == 0,  
            b S II - v II == 0, v II - g R == 0}, {S, II, R}]

\(\(sol\  = \ 
        NDSolve[{ode1, ode2, ode3, S[0] == 100, 
            II[0] == 1. , \ \ R[0] == 0}, {S[t], II[t], 
            R[t]}, {t, 0, 1000}];\)
\)
\(p1 = 
      Plot[Evaluate[{S[t]} /. sol], \ {t, 0, 
          500}, 
		PlotStyle -> {RGBColor[1, 0, 
              0]}, 
		PlotLabel -> "Susceptible"]

\)
\(p2 = 
      Plot[Evaluate[{II[t]} /. sol], \ {t, 0, 
          500}, 
		PlotStyle -> {RGBColor[1, 0, 
              0]}, 
		PlotLabel -> "Infected"]
\)
\(\(\(p3\)\(=\)\(Plot[
        Evaluate[{R[t]} /. sol], \ {t, 0, 
          500}, 
		PlotStyle -> {RGBColor[0, 0, 
              1]}, 
		PlotLabel -> "Recovered"]\)\(
\)\)\)

q1 = S[t] /. sol[[1]]
q2 = II[t] /. sol[[1]]
phaseplot2 = ParametricPlot[{q1, q2}, {t, 0, 500}]