external_pot = gp.MilkyWayPotential()
w0_1 = gd.PhaseSpacePosition(pos=[10, 0, 0] * u.kpc,
                             vel=[0, 200, 0] * u.km/u.s)
w0_2 = gd.PhaseSpacePosition(pos=w0_1.xyz + [10., 0, 0] * u.pc,
                             vel=w0_1.v_xyz + [0, 5, 0] * u.km/u.s)
w0 = gd.combine((w0_1, w0_2))
pot1 = gp.HernquistPotential(m=1e7*u.Msun, c=0.5*u.kpc, units=galactic)
particle_pot = [pot1, None]
nbody = DirectNBody(w0, particle_pot, external_potential=external_pot)
orbits = nbody.integrate_orbit(dt=1e-2*u.Myr, t1=0, t2=1*u.Gyr)

fig, ax = plt.subplots(1, 1, figsize=(5, 5)) # doctest: +SKIP
_ = orbits[:, 0].plot(['x', 'y'], axes=[ax]) # doctest: +SKIP
_ = orbits[:, 1].plot(['x', 'y'], axes=[ax]) # doctest: +SKIP