How to evaluate this integral analytically or numerically:
$$ I=\iiint dr_{12}dr_{13}dr_{14} $$
constrained by
$$ r_{12}<r_0,\\r_{23}<r_0,\\r_{34}<r_0, $$
where $r_0$ is a given real number and $$r_{ij}=\sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2}$$
I have been that this integral cannot be calculated analytically and one way for the numerical calculation is to use Monte Carlo methods, but I don't know how to implement this in MATLAB.