import gala.dynamics as gd

knots = np.linspace(0, 2, 21) * u.Gyr
masses = np.linspace(1e12, 4e12, 21) * u.Msun
pot = gp.TimeInterpolatedPotential(
    gp.HernquistPotential,
    time_knots=knots,
    m=masses,
    c=10.0 * u.kpc,
    units=galactic
)
w0 = gd.PhaseSpacePosition(
    pos=[8., 0, 0] * u.kpc,
    vel=[0, 220, 0] * u.km/u.s
)
orbit = gp.Hamiltonian(pot).integrate_orbit(
    w0, dt=1*u.Myr, n_steps=2000
)
fig = orbit.cylindrical.plot(["t", "rho"])
fig.axes[0].set_ylabel("radius $R$ [kpc]")