Appendix A of Liu, Baoyin, and Ma (2011) Equilibria, periodic orbits around equilibria, and heteroclinic connections in the gravity field of a rotating homogeneous cube shows an analytic expression for the gravitational potential of a uniform cube. I'd like to reproduce the orbit calculation using python, which will require evaluating the gradient of the reduced potential, i.e. acceleration.
Given sufficient coffee I can probably type that into Wolfram alpha, evaluate the gradient, and script as python. Alternately I could use a numerical gradient which would require four potential evaluations and be faster and easier for me to reliably script. I just want to make plots of single orbits, so I do not require very high accuracy.
Is there another way to get a fairly good approximation to the gravitational force from a uniform cube of side $a$ with similar (or less) effort? Fairly good might be say 1E-06 error at distances > $0.1a$ from a face, possibly need to stay farther from the corner for similar error, but it seems that most of the stable orbits tend to do that anyway.
3D direct integration at each time step is extremely slow with scipy's 'triple quadratic' method, not that I'd ever admit trying it. There may be an amazingly clever integration algorithm, and it could compete in speed considering I will be writing python without any of its numerical acceleration options at this point in time.