I managed to solve it using the following function:
given a cartesian point A and point B.
the geodesic path on a sphere is defined as:
r(t) = sin(1-t)*A + sin(t)*B, for t=[0, 1]
then normalize r(t)/||r(t)|| to get the points on a sphere
Can someone explain me why does it work?


What you have is a monotonic locus of $r(t)$ in both cases, however they are unrelated to the surface of the sphere. By normalising $r(t)$ afterwards to $\hat{r}(t)=\frac{r(t)}{\|r(t)\|}$, you're radially projecting your locus onto the unit sphere which gives you a non-linear parametrisation (non-linear interpolation). This non-linearity of parametrisation is visible in your plot. When you remove the sines there will be higher densities of points at the beginning and end of your locus than in the middle.