I would like to generate a random axis or unit vector in 3D. In 2D it would be easy, I could just pick an angle between 0 and 2*Pi and use the unit vector pointing in that direction. 

But in 3D I don't know how can I pick a random point on a surface of a sphere.

If I pick two angles the distribution won't be uniform on the surface of the sphere. There would be more points at the poles and less points at the equator.

If I pick a random point in the (-1,-1,-1):(1,1,1) cube and normalise it, then there would be more chance that a point gets choosen along the diagonals than from the center of the sides. So thats not good either.

But then what's the good solution? 

You need to use an equal-area projection of the sphere onto a rectangle.  Such projections are widely used in cartography to draw maps of the earth that represent areas accurately.

One of the simplest such projections is the axial projection of a sphere onto the lateral surface of a cylinder, as illustrated in the following figure:



This projection is area-preserving, and was used by Archimedes to compute the surface area of a sphere.

The result is that you can pick a random point on the surface of a unit sphere using the following algorithm:


Choose a random value of $\theta$ between $0$ and $2\pi$.
Choose a random value of $z$ between $-1$ and $1$.
Compute the resulting point:
$$
(x,y,z) \;=\; \left(\sqrt{1-z^2}\cos \theta,\; \sqrt{1-z^2}\sin \theta,\; z\right)
$$

