刚注意到置顶，请勿删。就是问个把图作出来的方法
p:=2.152 - 0.608 y[t] - 1.824 y[t]^2
q:=5.848 + 61.824 y[t] - 5.472 y[t]^2
NDSolve[{y'[t] + 3 z[t]*(y[t] + 1/6) == 0,
z'[t] + q/p*z[t]^2/(6 y[t] - 1) -
2 (1.152 - 1)*y[t] (6 y[t] - 1)/p + 50/p*(0.001)*(y[t] - 1/2) ==
0, y[0] == -0.1, z[0] == 0.1}, {y[t], z[t]}, {t, 0, 200}]

p := 2.152 - 0.608 y[t] - 1.824 y[t]^2
q := 5.848 + 61.824 y[t] - 5.472 y[t]^2
a = NDSolve[{y'[t] + 3 z[t]*(y[t] + 1/6) == 0,
z'[t] + q/p*z[t]^2/(6 y[t] - 1) -
2 (1.152 - 1)*y[t] (6 y[t] - 1)/p + 50/p*(0.001)*(y[t] - 1/2) ==
0, y[0] == -0.1, z[0] == 0.1}, {y[t], z[t]}, {t, 0, 200}]
PolarPlot[{y[t], z[t]} /. Flatten@a, {t, 0, 200}]