Conservation laws
I would start from the most basic question. Any operator evolves in time as:
$$O(t) = e^{\frac{i}{\hbar} Ht} O(0) e^{-\frac{i}{\hbar} Ht}$$
(this is a consequence of the Hamiltonian being the generator of time translations)
where $U(t) = e^{-\frac{i}{\hbar} Ht}$ is the time evolution operator. It follows straightforwardly that when we derive this equation with respect to time, we get:
$$\frac{d}{dt}O(t) = -\frac{i}{\hbar}\left[O(t),H\right]$$
therefore, if:
$$\left[O(t),H\right] = 0\rightarrow \frac{d}{dt}O(t) = 0$$
and so you can say that $O$ is conserved, i.e. it doesn't change with respect to time evolution.
Conservation of angular momentum
The conservation of angular momentum can be understood in terms of invariance under rotations (see edit below). While this explanation is quite elegant, it requires alot of prior knowledge, a more direct approach can be taken by considering the explicit definition of the orbital angular momentum $L$:
$$L_i = (x\times p)_i =\sum_{jk}\epsilon_{ijk}x_j p_k =-i\hbar\sum_{jk}\epsilon_{ijk}x_j\frac{\partial}{\partial x_k} $$
The Hamiltonian for a central potential is written as:
$$H = \frac{p^2}{2\mu}+U(r) = \frac{p^2_r}{2\mu}+\frac{L^2}{2\mu r^2}+U(r)$$
where $r = \sqrt{x^2+y^2+z^2}$.
Now, remembering that $[x_i,p_j] = i\hbar \delta_{ij}$ we have:
\begin{align}
[L_i,p^2] &= \sum_l [L_i,p_l]p_l+\sum_l p_l[L_i,p_l]\\
&=\sum_{ljk}\epsilon_{ijk}([x_j,p_l]p_kp_l+p_lp_k[x_i,p_l])\\
&=i\hbar \sum_{jk}\epsilon_{ijk}(p_jp_k+p_kp_j) = 0
\end{align}
since $\epsilon_{ijk}$ is antisymmetric.
Now, we use the know result $[p,f(x)]=-i\hbar\frac{\partial}{\partial x}f(x)$ to obtain:
\begin{align}
[L_i,U(r)] = \sum_{jk}\epsilon_{ijk}x_j[p_k,U(r)] =-i\hbar\sum_{jk}\epsilon_{ijk}x_j \frac{\partial r}{\partial x_k} \frac{\partial}{\partial r}U(r)
\end{align}
now:
$$\frac{\partial r}{\partial x_k} = \frac{x_k}{r} $$
therefore:
\begin{align}
[L_i,U(r)] = \sum_{jk}\epsilon_{jk}x_j[p_k,U(r)] =-i\hbar\frac{1}{r}\frac{\partial}{\partial r}U(r)\sum_{jk}\epsilon_{ijk}x_j x_k = 0
\end{align}
notice that $x_j x_k= \frac{1}{2}\left([x_j,x_k]+\{x_j,x_k\}\right) =\frac{1}{2}\{x_j,x_k\}$ since the coordinates commute among themeslves.
We conclude that: $$[L_i,H] = 0$$ and that: $$[L^2,H] = 0,\qquad [L_i,L^2] = 0$$
so the angular momentum is a conserved quantity. (As was specified in a comment this follows trivially from the rotational invariance of the Hamiltonian)
Factorization of eigenvalues
Commutatitvity implies that there exists a common basis of eigenvalues between $H$, $L_i$ and $L^2$ (commonly we choose $i=3$, i.e. $L_z$), we call these $\rvert n,l,m\rangle$:
\begin{align}
&H \rvert n,l,m\rangle= E_n \rvert n,l,m\rangle\\
&L_z \rvert n,l,m\rangle =\hbar m \rvert n,l,m\rangle\\
&L^2 \rvert n,l,m\rangle= \hbar^2 l(l+1)\rvert n,l,m\rangle
\end{align}
we define the wave functions in spherical coordinates:
\begin{align}
&x = r \sin\theta\cos\phi\\
&y =r \sin\theta\sin\phi\\
&z = r\cos\theta
\end{align}
as:
$$\psi_{nlm} = \langle r,\theta,\phi \rvert n,l,m\rangle=\psi_{nlm}(r,\theta,\phi)$$
now we project the first eigen-equation on the coordinate basis:
$$H\psi_{nlm}(r,\theta,\phi) =E_n\psi_{nlm}(r,\theta,\phi) = (\frac{p^2_r}{2\mu}+\frac{\hbar^2 l(l+1)}{2\mu r^2}+U(r))\psi_{nlm}(r,\theta,\phi)$$
now its obvious, this equation does not depend anymore on $\theta$ or $\phi$ and therefore the solution must simplify as:
$$\psi_{nlm}(r,\theta,\phi) = u_{nl}(r)Y_{lm}(\theta,\phi)$$
where $Y_{lm}(\theta,\phi)$ are the common eigenvalues of $L^2$ and $L_z$ just do the same on the eigen-equations of $L_i$ and $L^2$ and you'll see that:
\begin{align}
&L_z \psi_{nlm}(r,\theta,\phi) =\hbar m \psi_{nlm}(r,\theta,\phi)\\
&L^2 \psi_{nlm}(r,\theta,\phi)= \hbar^2 l(l+1)\psi_{nlm}(r,\theta,\phi)
\end{align}
where in the coordinate basis:
\begin{align}
&L_z = -i\hbar\frac{\partial}{\partial\phi}\\
&L^2 = -\hbar^2\left(\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\phi^2}\right)
\end{align}
and therefore variables separate in these equations too (only $Y_{lm}(\theta,\phi)$ remains).
We conclude that:
$$E_n u_{nl}(r) = (\frac{p^2_r}{2\mu}+\frac{\hbar^2 l(l+1)}{2\mu r^2}+U(r))u_{nl}(r)$$
Rotational invariance
Consider now a generic rotation of the coordinate axis representated by a matrix $R_{ij}$ such that:
$$x'_i = \sum_{j}R_{ij} x_j$$
and
$$p'_i = \sum_{j}R_{ij} p_j$$
it is clear that the square of any vector is left invariant under rotations since a rotation by definition cannot change its length:
$$x'^2 = \sum_i\sum_{j}R_{ij} x_j \sum_{k}R_{ik} x_k = x^2$$
then we find the condition:
$$\sum_{j} R_{ij}R_{jk} = \delta_{ik} $$
namely:
$$R^T R=1$$
where $R^T$ is the transpose matrix.
A rotation $R$ induces on the Hilbert space a unitary transformation $U(R)$ such that the coordinate operators tranform as:
$$U(R)^{\dagger}\hat{x}_i U(R) =\sum_j R_{ij}\hat{x}_j$$
this implies that our Hamiltonian is rotationally invariant:
$$U^{\dagger}(R)\,\hat{H}\, U(R) =\hat{H}$$
since it is only a function of $r$ and $p^2$.
Consider now a rotation by an infitesimal angle:
$R_{ij} = \delta_{ij}+\omega_{ij}+O(\omega^2)$
then it can be proven that:
$$U(1+\omega) = 1+ \frac{i}{2\hbar}\sum_{ij} \omega_{ij}L_{ij}$$
namely, the angular momentum is the generator of rotations. If we substitute this infinitsimal transformation into the tranformation law for the Hamiltonian we find:
$$[H,L] = 0$$
so the invariance of Hamiltonian under roations causes the angular momentum to be conserved.
