The method described in the article does not produce a spherically symmetric distribution since each value of z is equally likely. This means that a small circle near the north pole will contain as many points as the equator. Instead the different values for z should be selected with a probability that scales like 1/r where r is the radius of that circle (i.e. r = sin(theta)).

Alternatively you can generate random x, y, and z coordinates between -1 and 1 and toss it if it falls outside of the sphere (i.e. if x^2+y^2+z^2 > 1). Then just renormalized the coordinates so that they fall on the sphere.

The benefit of this is that the method easily generalizes to arbitrary dimension.