Convert Rotation matrix to Euler angles $~zyz~ (y$ convention$)$ analytically. The rotation matrix of Euler angle $ZYZ$ is:
$$
R_{z1}=\left[\cos(\psi),\sin(\psi),0;-\sin(\psi),\cos(\psi),0;0,0,1 \right];
$$
$$
R_y=[\cos(\theta),0,-\sin(\theta);0,1,0;\sin(\theta),0,\cos(\theta)];
$$
$$
R_{z2}=[\cos(\phi),\sin(\phi),0;-\sin(\phi),\cos(\phi),0;0,0,1];
$$
$R=Rz1*Ry*Rz2~$;
$$
R=\left[
\begin{array}{ccc}
\cos{\phi}\cos{\theta}\cos{\psi}-\sin{\phi}\sin{\psi} &\sin{\phi}\cos{\theta}\cos{\psi}+\cos{\phi}\sin{\psi}& -\sin{\theta} \cos{\psi}\\-\cos{\phi}\cos{\theta}\sin{\psi}-\sin{\phi}\cos{\psi}&-\sin{\phi}\cos{\theta}\sin{\psi}+\cos{\phi}\cos{\psi}& \sin{\theta}\sin{\psi}\\\cos{\phi}\sin{\theta}&\sin{\phi}\sin{\theta}&\cos{\theta}
\end{array}
\right]
$$
But if I have a rotation matrix R2 and I want to get Euler angle in $ZYZ~
(Y$ convention$)$ I can't find how to do it exactly.
for other Euler angles there is a solution, for example, see:
http://www.gregslabaugh.net/publications/euler.pdf
Why the following give me wrong answer when calculated in matlab:
$$ϕ=\text{atan}\left(\frac{R(3,2)}{R(3,1)}\right); ~~θ=\text{acos}(R(3,3)); ~~Ψ=\text{atan}\left(-\frac{R(2,3)}{R(1,3)}\right)$$
 A: The simplest approach to extract correctly Euler angles from a rotation matrix for any sequence of angles is using the $\mathrm{atan2}$ function. In the end, it is done in the same way (and maybe also explained why) in the text you linked. Note that compared to other inverse trigonometric functions, $\mathrm{atan2}$ has the range $(-\pi,\pi]$ (the full circle). 
In your case of the $z$-$y$-$z$ rotation, if $\sin\theta\neq 0$, then
$$
\phi=\mathrm{atan2}(R_{32},R_{31}), \quad \psi=\mathrm{atan2}(R_{23},-R_{13}).
$$
We can obtain the $\theta$ angle from the last row or the last column of $R$. For example, considering the last row of $R$, we have
$$
R_{31}\cos\phi+R_{32}\sin\phi=\sin\theta(\cos^2\phi+\sin^2\phi)=\sin\theta
$$
so
$$
\theta=\mathrm{atan2}(R_{31}\cos\phi+R_{32}\sin\phi,R_{33}).
$$
For these Euler angles to be well-defined, the condition $\sin\theta\neq 0$ is required (that is, $\theta\neq k\pi$, $k\in\mathbb{Z}$). Otherwise, e.g., if $\theta=0$, the two $z$-angles are not uniquely defined.
A: When you know the tangent of an angle, there is an indeterminacy of a half-turn on the angle.
I recommend to avoid the use of the arc cosine, and prefer the arc tangent, with
$$\tan\theta=\frac{\sqrt{(\cos\phi\sin\theta)^2+{(\sin\phi\sin\theta)^2}}}{\cos\theta}.$$ Beware that this drops the sign of $\sin\theta$.
After obtaining the angles, plug them back to compare to the original  matrix and adjust for the quadrants. (Sorry, I don't have time to deal with the complete discussion.)
