Next, we will calculate the commutation relations between the different angular momentum operators, to determine which of them can be known simultaneously.
\begin{equation} [\hat{l}_x , \hat{l}_y]=\hat{l}_x \hat{l}_y - \hat{l}_y \hat{l}_x \end{equation}
So that expressions are not too long, let's calculate the two terms separately
\begin{equation} \hat{l}_x\hat{l}_y f=\left[-i\hbar\left(y\frac{\partial} {\partial z}-z\frac{\partial}{\partial y}\right)\right]\left[-i\hbar\left(z\frac{\partial f}{\partial x}-x\frac{\partial f}{\partial z}\right)\right]= \end{equation}
\begin{equation} =-\hbar^2\left(y\frac{\partial f}{\partial x} +yz\frac{\partial^2 f}{\partial z \partial x}-yz\frac{\partial^2 f}{\partial z^2}-z^2\frac{\partial^2 f} {\partial y \partial x}+zx\frac{\partial^2 f}{\partial y \partial z}\right) \end{equation}
Now we calculate the second term of the commutator:
\begin{equation} \hat{l}_y\hat{l}_x f=\left[-i\hbar\left(z\frac{\partial}{\partial x}-x\frac{\partial}{\partial z}\right)\right]\left[-i\hbar\left(y\frac{\partial f}{\partial z}-z\frac{\partial f}{\partial y}\right)\right]=\end{equation}
\begin{equation} =-\hbar^ 2 \left( zy\frac{\partial^2 f}{\partial x \partial z}-z^2\frac{\partial^2 f}{\partial x \partial y}-zy\frac{\partial ^2 f}{\partial z^2}+x\frac{\partial f}{\partial y}+xz\frac{\partial^2 f}{\partial z \partial y}\right) \end{equation}
Substituting in the commutator:
\begin{equation} [\hat{l}_x , \hat{l}_y]=\hat{l}_x \hat{l}_y - \hat{l}_y \hat{ l}_x =-\hbar^2\left(y\frac{\partial}{\partial x}-x\frac{\partial}{\partial y}\right)=-\frac{\hbar}{i }\left(-i\hbar\right)\left(x\frac{\partial}{\partial y}-y\frac{\partial}{\partial x}\right)=i\hbar\hat{l }_z \end{equation}
The commutators between the remaining angular momentum components are calculated analogously. Summarizing:
\begin{eqnarray} \left[\hat{l}_x , \hat{l}_y\right] &=& i\hbar\hat{l}_z\\ \left[\hat{l}_y , \hat{l}_z\right] &=& i\hbar\hat{l}_x\\ \left[\hat{l}_z , \hat{l}_x\right] &=& i\hbar\hat {l}_y \end{eqnarray}
Let's evaluate the commutators of $\hat{l}^2$
\begin{equation} \left[\hat{l}^2, \hat{l}_x \right]=\left[\hat{l}^2_x + \hat{l}^2_y + \hat{l} ^2_z ,\hat{l}_x\right]=\underbrace{\left[\hat{l}^2_x ,\hat{l}_x\right]}_{0}+\left[\hat{l} ^2_y , \hat{l}_x\right]+\left[\hat{l}^2_z , \hat{l}_x\right]=\left[\hat{l}_y\hat{l}_y , \hat{l}_x\right]+\left[\hat{l}_z\hat{l}_z , \hat{l}_x\right]= \end{equation}
Applying the property $\left[\hat{A}\hat{B},\hat{C}\right]=\left[\hat{A},\hat{C}\right]\hat{B} + \hat{A}\left[\hat{B},\hat{C}\right]$
\begin{equation} =\left[\hat{l}_y,\hat{l}_x \right]\hat{l}_y + \hat{l}_y\left[\hat{l}_y,\hat {l}_x \right]\hat{l}_y + \left[\hat{l}_z,\hat{l}_x \right]\hat{l}_z + \hat{l}_z\left[\ hat{l}_z,\hat{l}_x \right]\hat{l}_z =\\-i\hbar\hat{l}_z\hat{l}_y + \hat{l}_y (-i \hbar\hat{l}_z) +i\hbar\hat{l}_y \hat{l}_z + \hat{l}_z(i\hbar\hat{l}_y)=0 \end{equation}