epsilon = 0.1
tMin = 0
tMax = 20
T.ode_solve(t_span = [tMin, tMax], params=(epsilon, ), num_points=1000)
(tN, YN) = zip(*T.solution)
(yN, ypN) = zip(*YN)
PN = line(zip(tN, yN), color = 'blue')
PR0 = plot(yR0, (t, tMin, tMax), color = 'red')
PR1 = plot(yR0 + epsilon * yR1, (t, tMin, tMax), color = 'green')
PTT0a = plot(yTT0a(epsilon = epsilon), (t, tMin, tMax), color = 'orange')
PTT0ab = plot(yTT0ab(epsilon = epsilon), (t, tMin, tMax), color = 'purple')
show(PN + PR0 + PR1 + PTT0a + PTT0ab)