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. When $n\rightarrow\infty$, the field should approach the real distribution.

I am rewriting my comment as an answer.

I think that you have several options. Some of them are:

• 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 Chebyshev 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. When $n\rightarrow\infty$, the field should approach the real distribution.

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.

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. When $n\rightarrow\infty$, the field should approach the real distribution.