dxyz = orbits[:, 0].xyz - orbits[:, 1].xyz

fig, ax = plt.subplots(1, 1, figsize=(5, 5)) # doctest: +SKIP
ax.plot(dxyz[0], dxyz[1]) # doctest: +SKIP
ax.set_xlim(-0.1, 0.1)
ax.set_ylim(-0.1, 0.1)