I'm rewriting my comment as an answer.

I think that you have some options:
 
 - You can use [SymPy]( http://docs.sympy.org/latest/tutorial/index.html) to find the gravitational field from the potential. Then you can generate your [Python code from it]( http://docs.sympy.org/latest/modules/printing.html).

 
 - 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]( https://github.com/pychebfun/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. When $n\rightarrow\infty$, the field should approach the real distribution.