def METHOD_NAME(fig, axes, K_values, theta_values):
    '''
    Makes animated subplots in the given ``axes``, where each ``theta_values`` entry
    is the full recording of ``Theta`` from the monitor.
    '''
    artists = []
    for ax, K, Theta in zip(axes, K_values, theta_values):
        x, y = to_x_y(Theta.T[0])
        dots = ax.scatter(x, y, c=Theta.T[0])
        r, phi = calc_coherence_and_phase(x, y)
        arrow = ax.arrow(0, 0, r*np.cos(phi), r*np.sin(phi), color='black')
        mean_dot, = ax.plot(r*np.cos(phi), r*np.sin(phi), 'o', color='black')
        if abs(K) > 0:
            title = f"coupling strength K={K:.1f}"
        else:
            title = "uncoupled"
        ax.text(-1., 1.05, title, color='gray', va='bottom')
        ax.set_aspect('equal')
        ax.set_axis_off()
        ax.set(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2))
        artists.append((dots, arrow, mean_dot))
    def update(frame_number):
        updated_artists = []
        for (dots, arrow, mean_dot), K, Theta in zip(artists, K_values, theta_values):
            t = frame_delay*frame_number
            ts = timestep(t, defaultclock.dt)
            x, y = to_x_y(Theta.T[ts])
            dots.set_offsets(np.vstack([x, y]).T)
            r, phi = calc_coherence_and_phase(x, y)
            arrow.set_data(dx=r*np.cos(phi), dy=r*np.sin(phi))
            mean_dot.set_data(r*np.cos(phi), r*np.sin(phi))
            updated_artists.extend([dots, arrow, mean_dot])
        return updated_artists
    ani = animation.FuncAnimation(fig, update, frames=int(magic_network.t/frame_delay),
                                interval=20, blit=True)
    return ani