You can also do this in spherical coordinates, in which case there is no rejection. First you generate the radius and the two angles at random, then you use the transition formula to recover $x$, $y$ and $z$ ($x = r \sin \theta \cos \phi$, $y = r \sin \theta \sin \phi$, $z = r \cos \theta$).
You generate $\phi$ unifomly between $0$ and $2\pi$. The radius $r$ and the inclination $\theta$ are not uniform though. The probability that a point is inside the ball of radius $r$ is $r^3$ so the probability density function of $r$ is $3 r^2$. You can easily check that the cubic root of a uniform variable has exactly the same distribution, so this is how you can generate $r$. The probability that a point lies within a spherical cone defined by inclination $\theta$ is $(1-\cos\theta)/2$ or $1 - (1-\cos (-\theta))/2$ if $\theta > \pi/2$. So the density $\theta$ is $sin(\theta)/2$. You can check that minus the arccosine of a uniform variable has the proper density.
Or more simply, we can simulate the cosine of $\theta$ uniformly beteen $-1$ and $1$.
In R this would look as shown below.
n <- 10000 # For example n = 10,000.
phi <- runif(n, max=2*pi)
r <- runif(n)^(1/3)
cos_theta <- runif(n, min=-1, max=1)
x <- r * sqrt(1-cos_theta^2) * cos(phi)
y <- r * sqrt(1-cos_theta^2) * sin(phi)
z <- r * cos_theta
In the course of writing and editing this answer, I realized that the solution is less trivial that I thought.
I think that the easiest and computationally most efficient method is to follow @whuber's method to generate $(x,y,z)$ on the unit sphere as shown on this post and scale them with $r$.
xyz <- matrix(rnorm(3*n), ncol=3)
lambda <- runif(n)^(1/3) / sqrt(rowSums(xyz^2))
xyz <- xyz*lambda