pot = gp.MiyamotoNagaiPotential(m=1E11, a=6.5, b=0.27, units=galactic)

fig, ax = plt.subplots(1,1) # doctest: +SKIP
pot.plot_contours(grid=(np.linspace(-15,15,100), 0., 1.), marker='', ax=ax) # doctest: +SKIP
E_unit = pot.units['energy'] / pot.units['mass']
ax.set_xlabel("$x$ [{}]".format(pot.units['length'].to_string(format='latex'))) # doctest: +SKIP
ax.set_ylabel("$\Phi(x,0,1)$ [{}]".format(E_unit.to_string(format='latex'))) # doctest: +SKIP
fig.tight_layout()