import astropy.units as u
import gala.dynamics as gd
import gala.potential as gp
from gala.units import galactic

pot = gp.PlummerPotential(m=1E10 * u.Msun, b=1. * u.kpc, units=galactic)
w0 = gd.PhaseSpacePosition(pos=[2.,0,0] * u.kpc,
                           vel=[0.,75,15] * u.km/u.s)
orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=1., n_steps=5000)
fig = orbit.plot()