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.
