Your missing the requirement that the random sample be uniform in the sphere. That is, if the sphere is cut into any number of pieces of equal volume, the sample is equally likely to appear in any piece.
Your approach samples uniformly from the radius, so you're going to end up with far too many samples close to the center.
Interesting... thanks for the pointer. I did some testing and the points definitely do seem to cluster along the axis... just increasing the number of iterations of √(r²-n²) seemed to improve it though. At 5, I can't tell the difference anymore:
Your approach samples uniformly from the radius, so you're going to end up with far too many samples close to the center.