epsilon = 0.3
tMax = 10
# Numerically solve
solver.ode_solve(params = [epsilon], t_span = [0, tMax], num_points = 1000)
# Plot u
(plot(u0, (t, 0, tMax), color = 'red')
+ plot(u0 + epsilon * u1, (t, 0, tMax), color = 'green')
+ line(((s[0], s[1][0]) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'u'))
# Plot 1
(plot(q0, (t, 0, tMax), color = 'red')
+ plot(q0 + epsilon * q1, (t, 0, tMax), color = 'green')
+ line(((s[0], s[1][1]) for s in solver.solution), color = 'blue')).show(axes_labels = ('t', 'q'))