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
    wdot = np.zeros_like(w)
    wdot[0] = p
    wdot[1] = -np.sin(q) + A*omega_D*np.cos(omega_D*t)
    return wdot

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)))