Evaluating gradient of numerically approximated function

edit

Hello, I'd like to ask how to solve the following problem in Sage.

Let's suppose that I have a scalar function on R^3 (written in spherical coordinates) V(r,\theta,\phi), given by some tedious integral over a region in R^2: V(r,th,ph) = \int dx \int dy f(x,y,r,th,ph) (typical when calculating potentials by Green functions). The first question, how could I define such a function V in Sage? I guess it would involve approximating this function point by point by a 3-dim matrix containing results of numerical integration V(r0,th0,ph0) at different (r0,th0,ph0)...

Then, I would like to calculate the gradient of this V(r,th,ph) and plot streamlines or at least some plot of resulting vector field. I know how to do this when Sage is able to evaluate this gradient symbolically, but I don't know how to approximate it numerically and then make a plot.
I appreciate any help.