I'm rewriting my comment as an answer.
I think that you have some options:
You can use SymPy to find the gravitational field from the potential. Then you can generate your Python code from it.
You could create some kind of mesh and then compute the gradient numerically using finite differences, piecewise polynomials (FEM-like), or Chebysev polynomials (see pychebfun).
You can use a Boundary Integral Representation for the gravitational field of your problem.
You can represent your cube as a set of $n$ discrete particles and consider that the overall field is just the superposition of these "little" contributions. Even $n\rightarrow\infty$, the field should approach the real distribution.