You should be able to get accurate results with mpmath, a Python module for arbitrary-precision floating-point computations. There are examples of integration with singularities in the documentation. You'll want to explicitly tell it to break up the interval:

from mpmath import *
f = lambda x,y,z: 1./(x**2+y**2+z**2)**1./3
quad(f,[-1,0,1],[-1,0,1],[-1,0,1])

You may need to increase the precision (e.g. mp.dps=30) and it will likely be slow, but should be quite accurate.

You could also try nesting calls to MATLAB's quadgk, which uses adaptive Gauss-Kronrond quadrature in 1D.

You should be able to get accurate results with mpmath, a Python module for arbitrary-precision floating-point computations. There are examples of integration with singularities in the documentation. You'll want to explicitly tell it to break up the interval:

from mpmath import *
f = lambda x,y,z: 1./(x**2+y**2+z**2)**1./3
quad(f,[-1,0,1],[-1,0,1],[-1,0,1])

You may need to increase the precision (e.g. mp.dps=30) and it will likely be slow, but should be quite accurate.

You could also try nesting calls to MATLAB's quadgk(), which uses adaptive Gauss-Kronrod quadrature in 1D.