Comment (too long for comment format): Your outline toward an analytic solution seems feasible. I haven't tried it, but it might be easier to stay with rectangular coordinates, initially integrating over a uniform distribution on the ball.
Below are some simulated distributions, means, and standard deviations.
x1 is one coordinate;
R is squared distance from origin;
r is distance from origin.
m = 10^6
u1 = runif(m,-1,1); u2 = runif(m,-1,1); u3 = runif(m,-1,1)
ball = (u1^2 + u2^2 + u3^2 <= 1)
 523552 # points inside ball
 0.5235988 # volume of ball
# eliminate points outside ball
x1 = u1[ball]; x2 = u2[ball]; x3 = u3[ball]
 -0.0001425727 # aprx E(X1) = 0
hist(x1, prob=T, col="skyblue2")
# for clarity: reduced nr of pts in scatterplot
R = x1^2 + x2^2 + x3^2
hist(R, prob=T, col="skyblue2")
r = sqrt(R)
hist(r, prob=T, col="skyblue2")