Next, we will calculate the commutation relations between the different angular momentum operators, to determine which of them can be known simultaneously.

$$[\hat{l}_x , \hat{l}_y]=\hat{l}_x \hat{l}_y - \hat{l}_y \hat{l}_x$$

So that expressions are not too long, let's calculate the two terms separately

$$\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]=$$

$$=-\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)$$

Now we calculate the second term of the commutator:

$$\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]=$$

$$=-\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)$$

Substituting in the commutator:

$$[\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$$

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$

$$\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]=$$

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]$

$$=\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$$