Inspired by this recent Numberphile video, here is a demonstration of chaos in a simple dynamical system: two balls, with near-identical starting conditions, bounce around elastically off a circular wall. After a short time, the balls' trajectories diverge completely.
solve_ivp can be used to integrate the relevant equations of motion. Each part of the trajectory between bounces is integrated in turn, with the impact with the wall identified with an event occuring when the ball's position is greater than the radius of the circle. When the event is triggered, the ball's velocity is reflected (in the line of the normal to the wall at the impact point) and the next integration started.
import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt import matplotlib.animation as animation # Acceleration due to gravity, m.s-2 (downwards!). g = 10 # Circle radius, m. R = 1 # Time step, s. dt = 0.001 def solve(u0): """Solve the equation of motion for a ball bouncing in a circle. u0 = [x0, vx0, y0, vy0] are the initial conditions (position and velocity). """ # Initial time, final time, s. t0, tf = 0, 10 def fun(t, u): """Return the derivatives of the dynamics variables packed into u.""" x, xdot, y, ydot = u xddot = 0 yddot = -g return xdot, xddot, ydot, yddot def event(t, u): """If the ball hits the wall of the circle, trigger the event.""" return np.hypot(u, u) - R*1.01 # Make sure the event terminates the integration. event.terminal = True # Keep track of the ball's position in these lists. x, y = ,  while True: # Solve the equations until the ball hits the circular wall or until # the time tf. soln = solve_ivp(fun, (t0, tf), u0, events=event, dense_output=True) if soln.status == 1: # We hit the wall: save the path so far... tend = soln.t_events nt = int(tend - t0 / dt) + 1 tgrid = np.linspace(t0, tend, 100) sol = soln.sol(tgrid) x.append(sol) y.append(sol) # ...and restart the integration with the reflected velocities as # the initial conditions. u = soln.y[:, -1].copy() p = np.array((u, u)) p = p / np.linalg.norm(p) v = np.array((u, u)) v = v - 2 * (v @ p) * p u0 = p, v, p, v t0 = soln.t[-1] else: # We're done up to tf (or, rather, the last bounce before tf). break # Concatenate all the paths between bounces together. return np.concatenate(x), np.concatenate(y) # For the animation, set up the path lines and circle patches representing the # balls. fig, ax = plt.subplots() ax.set_xlim(-1, 1) ax.set_ylim(-1, 1) # Make the circular wall a little bit bigger because the ball paths eat into it. ax.add_patch(plt.Circle((0, 0), 1.02, fc='k', ec='w')) ax.axis('equal') ax.axis('off') line0, = ax.plot(, ) line1, = ax.plot(, ) pos =  u0 = [0.001, 0, 0, 0] ball0 = ax.add_patch(plt.Circle((u0, u0), 0.05, fc='tab:blue', ec='none')) pos.append(solve(u0)) u0 = [0.0015, 0, 0, 0] ball1 = ax.add_patch(plt.Circle((u0, u0), 0.05, fc='tab:orange', ec='none')) pos.append(solve(u0)) def init(): """ Initialization, because we're blitting and need references to the animated objects. """ return line0, line1, ball0, ball1 def animate(i): """Draw frame i of the animation.""" line0.set_data(pos[:i], pos[:i]) ball0.set_center((pos[i], pos[i])) line1.set_data(pos[:i], pos[:i]) ball1.set_center((pos[i], pos[i])) return line0, line1, ball0, ball1 interval, nframes = 1000 * dt, int(len(pos)) ani = animation.FuncAnimation(fig, animate, frames=nframes, repeat=False, init_func=init, interval=interval, blit=True) plt.show()