import astropy.units as u
import gala.potential as gp
import gala.dynamics as gd
from gala.units import galactic
pot = gp.MiyamotoNagaiPotential(m=2.5E11, a=6.5, b=0.26, units=galactic)
w0 = gd.PhaseSpacePosition(pos=[11., 0., 0.2]*u.kpc,
                           vel=[0., 200, 100]*u.km/u.s)
orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=1., n_steps=1000)
fig = orbit.plot()