The Hamiltonian of spin-orbit coupling in an external magnetic field $\vec{B} = B\vec{e}_z$ is given by $$ H = \beta L\cdot S+\frac{\mu_b}{\hbar}(L_z+2S_z)B.$$ The ladder operators are $$ L_\pm = L_x \pm iL_y\qquad S_\pm = S_x \pm iS_y$$ and operate on the angular momentum eigenstates as follows:
$$ L_\pm|l,m_l\rangle = \hbar\sqrt{(l\pm m_l+1)(l\mp m_l)}|l,m_l\pm 1\rangle$$ $$S_\pm|s,m_s\rangle = \hbar\sqrt{(s\pm m_s+1)(s\mp m_s)}|s,m_s\pm 1\rangle$$
The Hamiltonian can then be expressed in terms of these operators:
$$ H = \frac \beta 4 (L_+S_-+L_-S_+)+L_zS_z+\frac{\mu_B}{\hbar} (L_z+2S_z)$$
My question is how to calculate the matrix elements of $H$ using the uncoupled basis states $|m_l,m_s\rangle$, so $\langle m_l',m_s'|H|m_l,m_s\rangle$, for a $p$-electron ($l=1$). I know the result which I found in a book, but that book doesn't show the explicit computation of the matrix entries. Since I'm quite new to this kind of things, it would be great if someone could show me how to do the calculation. I've been trying things for days now but nothing turned out to be correct.