I can't wrap my head around a problem, I have a methane molecule, which is basically a tetrahedron, and I orient one of its bonds (let's say from atom 1 and atom 2) along the $z$ axis of a cubic box.
Now I know that the probability distribution for the angle $A$ between the 1-2 vector and the $z$ axis would be shaped as a sine function, but i think I understood that $\cos(A)$ instead should have a flat probability $P[\cos(A)]$ distribution (from https://www.lockhaven.edu/~dsimanek/puzzles/hatbox.htm#:~:text=While%20puzzling%20how%20to%20find,radius%20R%20and%20height%202R.%22) .
Now If I make my molecule rotate around the $y$ axis, so that I have 360 structures all equispaced of 1 degree I expected to see a flat $P[\cos(A)]$ distribution, and instead i get a distribution shaped like this:
Then I tried to get random numbers between 0 and 360 and use them to rotate the molecule around $x, y$ and $z$ of that number of degrees, and I still get the same distribution.
Then I tried to get three random numbers between $-1$ and $1$ and using arccos of those numbers to rotate the molecule around $x,y$ and $z$. And I still get the same shape of the distribution.
Can someone help me wrap my head around this?

What you are looking for is the distribution of the random variable $Y=\cos (X)$, where the angle $X$ is a random variable with the probability distribution function $$f_X(x)=\frac1\pi\quad\text {for}\quad 0 \le x \le \pi. $$ [Note that we need not to consider the range $[\pi,2\pi]$ in view of $\cos(2\pi-x)=\cos x$.]
Let us solve the problem using the distribution function technique. First we find the cumulative probability distribution function:
$$F_Y (y)=P (Y\le y)=P (\cos (X)\le y)=P (X\ge \arccos (y))\\ =\int\limits_{\arccos(y)}^\pi f_X (x)dx=1-\frac1\pi\arccos (y) .$$ The probability density function is then: $$f_Y (y)=\frac d {dy}F_Y(y)=\frac1 {\pi\sqrt {1-y^2}}, $$ in (at least qualitative) agreement with the result of your numerical simulation.