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.0, 1.0), 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()