import astropy.units as u
import matplotlib.pyplot as pl
import numpy as np
import gala.integrate as gi

def F(t, w, A, omega_D):
    q,p = w
    q_dot = p
    p_dot = -np.sin(q) + A*omega_D*np.cos(omega_D*t)
    return np.array([q_dot, p_dot])

integrator = gi.DOPRI853Integrator(F, func_args=(0.07, 0.75))
orbit = integrator.run([3.,0.], dt=0.1, n_steps=10000)
fig = orbit.plot(subplots_kwargs=dict(figsize=(8,4)))