First, I think your $H_{SOC}$ is incorrect and should read:
$$
H_{SOC}=\frac{\alpha}{2}\left(L_+\sigma_- + L_-\sigma_+ + L_z\sigma_z\right)\, .
\tag{1}
$$
The twist with your problem is that it uses the real spherical harmonics, rather than the complex exponential form more common in physics and better adapted to the evaluation of matrix elements. The conversion between the two can be found on this wiki page.
Thus:
\begin{align}
p_y&=\frac{i}{\sqrt{2}}\left(Y_{1}^{-1}+Y_1^1\right)\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{2}\\
p_x&=\frac{1}{\sqrt{2}}\left(Y_{1}^{-1}-Y_1^1\right)
\quad\Rightarrow \quad \vert p_y\rangle = \frac{i}{\sqrt{2}}\left(\vert 1,-1\rangle + \vert 1,1\rangle\right)\, ,\tag{3}\\
p_z&=Y_1^0\qquad\qquad \qquad\quad\Rightarrow \quad \vert p_z\rangle = \vert 1,0\rangle\, , \tag{4}
\end{align}
where the arguments $(\theta,\varphi)$ of the $Y_{\ell}^{m}$ have been omitted for clarity.
The next step is to recall the action of the various operators on $Y_{\ell m}$:
\begin{align}
\hat L_\pm Y_{\ell}^{m}&=\sqrt{(\ell\mp m)(\ell\pm m+1)}Y_{\ell}^{m\pm 1}\, ,\tag{5}\\
\hat L_z Y_{\ell}^{m}&= mY_{\ell}^{m}\tag{6}
\end{align}
from which one obtains for instance (if I did not make obvious mistakes), the matrix representation of $\hat L_+$ in the basis $\{\vert p_x\rangle, \vert p_y\rangle, \vert p_z\rangle\}$
\begin{align}
\hat L_+\vert p_x\rangle &= \vert p_z\rangle \qquad
\hat L_+\vert p_y\rangle = i \vert p_z\rangle \qquad
\hat L_+\vert p_z\rangle =-\vert p_x\rangle - i\vert p_y\rangle\, , \tag{7}
\end{align}
so that
\begin{align}
\hat L_+\to \left(\begin{array}{ccc}
0&0&-1\\
0&0&-i\\
1&i&0\end{array}\right)\, ,\qquad \hat L_-\to
\left(\begin{array}{ccc}
0&0&1\\
0&0&-i\\
-1&i&0\end{array}\right)\, . \tag{8}
\end{align}
since the matrix representation for $\hat L_-$ will be transpose conjugate of the matrix representation for $\hat L_+$. Likewise:
\begin{align}
\hat \sigma_+\vert \uparrow \rangle &=0\, , \quad
\hat \sigma_+\vert \downarrow \rangle =\vert \uparrow\rangle\, ,\quad
\hat \sigma_-\vert \uparrow \rangle =\vert\downarrow\rangle\, , \quad
\hat \sigma_-\vert \downarrow \rangle =0\, . \tag{9}
\end{align}
Hence:
\begin{align}
\hat L_+\hat \sigma_- \vert p_x;\uparrow\rangle &= \Bigl[\hat L_+\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_-\vert \uparrow\rangle\Bigr]= \vert p_z\rangle \vert\downarrow\rangle= \vert p_z;\downarrow\rangle\, , \tag{10}\\
\hat L_-\hat \sigma_+ \vert p_x;\uparrow\rangle &= 0 \, ,\tag{11} \\
\hat L_z\hat\sigma_z\vert p_x;\uparrow\rangle &=
\Bigl[\hat L_z\vert p_x\rangle \Bigr]
\Bigl[\hat \sigma_z\vert \uparrow\rangle\Bigr]=i\vert p_y\rangle \vert \uparrow\rangle=
i\vert p_y;\uparrow\rangle \tag{12}
\end{align}
This is not quite the first column of your matrix, but as I understand it this matrix is not quite correct either since it ought to be hermitian but is not: for instance the matrix element in position $(2,6)$ should be the complex conjugate of the matrix element in position $(6,2)$.
Anyways, my calculation would give the elements in the first column as
$$
(0,0,i,0,0,1)^T \tag{13}
$$
This might be because the basis elements are ordered as
$\{\vert p_x\uparrow\rangle,\vert p_y\uparrow\rangle,\vert p_z\uparrow\rangle,\vert p_x\downarrow\rangle,\vert p_y\downarrow\rangle,\vert p_z\downarrow\rangle\}$ rather than the ordering you give (or else it's quite possible I made an error somewhere). With this ordering my first column would come out as $(0,i,0,0,0,1)^T$, which is identical to yours. Note that this does not fix the hermiticity problem mentioned earlier.
The other columns are found in the same manner.
In addition to the answer by hsinghal it is worth point out some historical notational quirks.
An expression such as $^2P_{3/2}$ is called a Term Symbol. The superscript is the multiplicity of the electron spins, i.e. 2$S$+1 for total spin S. The capital letter, P in this example is the total orbital angular momentum and has letters and values of S=0, P=1, D=2, etc (as for atomic orbitals) and the subscript is the total spin and orbital angular momentum given by J. Thus the term symbol is $^{2S+1}L_J$
It should also be noted that when calculating values such as $L= |l_1−l_2| $ to $|l_1+l_2|$ the values are taken in unit steps, i.e. each values differs by 1 from the next until no more can fit with the formula. Sometimes there is only one value.
The origin of the Zeeman effect is in the magnetic field produced by a 'spinning' charge. In the absence of an external field the sub-levels of J (2$J$+1 of them) have the same energy. In an external magnetic field the magnetic field caused by spin interacts with the external field and levels split. (The hamiltonian is $H=-\mu . B$, where $\mu$ is the magnetic dipole of the atom and is directly proportional to J, and B is the external magnetic field strength). You can use Hund's Rules to decide which state is highest or lowest. See chapter 2 in 'Molecules and Radiation' by J. Steinfeld. The figure shows an example of a $p^2$ configuration.
Nowadays the Zeeman effect is tremendously important, not in atomic spectroscopy but in nuclear magnetic resonance (NMR). This is widely used in Chemistry to determine the structure of molecules and even of proteins when in solution. NMR is the basis of MRI imaging widely used in hospitals. In this case it is the nuclear Zeeman effect, usually using spin 1/2 protons, but also very many other nuclei with integer and half integer spin can be used. See 'Spin Dynamics' M. Levitt.
Best Answer
The spin-orbit interaction is responsible for splitting the $5p$ level in two, hence the two identified $5s \rightarrow 5p$ wavelengths. To calculate the splitting magnitude, you want to be looking at the difference $\Delta E$ between the two transition energies.
By the way, I'm seeing the coupling constant typically characterized in inverse centimeters: $$ a = \frac{2}{3} \frac{\Delta E}{hc} $$