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.
First let's add the first two spins $\mathbf{s}_1+\mathbf{s}_2=\mathbf{s}_{12}$, this can occur in two ways: $s_{12}=0$ and $s_{12}=1$. The general rule of momentum coupling is
$$
|s_1s_2:s_{12}m_{12}\rangle_{12}=\sum_{m_1m_2}( s_1m_1s_2m_2|s_{12}m_{12})|s_1m_1\rangle_1|s_2m_2\rangle_2.
$$
In the case of two spin-one-half particles it leads to familiar singlet and triplet states:
$$
{\textstyle\left|\frac12\frac12:00\right\rangle_{12}}=\frac{|\uparrow\rangle_1|\downarrow\rangle_2-|\downarrow\rangle_1|\uparrow\rangle_2}{\sqrt2},
$$
$$
{\textstyle\left|\frac12\frac12:11\right\rangle_{12}}=|\uparrow\rangle_1|\uparrow\rangle_2,\quad
{\textstyle\left|\frac12\frac12:10\right\rangle_{12}}=\frac{|\uparrow\rangle_1|\downarrow\rangle_2+|\downarrow\rangle_1|\uparrow\rangle_2}{\sqrt2},\quad{\textstyle\left|\frac12\frac12:1-1\right\rangle_{12}}=|\downarrow\rangle_1|\downarrow\rangle_2
$$
(in the following we will omit $\frac12\frac12:$ in these state vectors).
Now let's add the third spin $s_3=\frac12$ to $s_{12}=0,1$, the resulting total spin is $s$ and its projection is $m$. It occurs according to the rule
$$
|s_1s_2(s_{12})s_3:sm\rangle=\sum_{m_{12}m_3}(s_{12}m_{12}s_3m_3|sm)|s_{12}m_{12}\rangle_{12}|s_3m_3\rangle_3.
$$
The possible cases are:
A) $s_{12}=0$, $s=\frac12$, $m=\pm\frac12$:
$$
\textstyle\left|\frac12\frac12(0)\frac12:\frac12\frac12\right\rangle=\left|00\right\rangle_{12}\left|\uparrow\right\rangle_3.
$$
$$
\textstyle\left|\frac12\frac12(0)\frac12:\frac12-\frac12\right\rangle=\left|00\right\rangle_{12}\left|\downarrow\right\rangle_3.
$$
B) $s_{12}=1$, $s=\frac12$, $m=\pm\frac12$:
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac12\frac12\right\rangle}=\sqrt{\frac23}{\textstyle\left|11\right\rangle_{12}\left|\downarrow\right\rangle_3}-\sqrt{\frac13}{\textstyle\left|10\right\rangle_{12}\left|\uparrow\right\rangle_3},
$$
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac12-\frac12\right\rangle}=
\sqrt{\frac13}{\textstyle\left|10\right\rangle_{12}\left|\downarrow\right\rangle_3}-\sqrt{\frac23}{\textstyle\left|1-1\right\rangle_{12}\left|\uparrow\right\rangle_3},
$$
C) $s_{12}=1$, $s=\frac32$, $m=\pm\frac12,\pm\frac32$:
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac32\frac32\right\rangle}=
{\textstyle\left|11\right\rangle_{12}\left|\uparrow\right\rangle_3},
$$
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac32\frac12\right\rangle}=
\sqrt{\frac13}{\textstyle\left|11\right\rangle_{12}\left|\downarrow\right\rangle_3}+\sqrt{\frac23}{\textstyle\left|10\right\rangle_{12}\left|\uparrow\right\rangle_3},
$$
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac32-\frac12\right\rangle}=
\sqrt{\frac23}{\textstyle\left|10\right\rangle_{12}\left|\downarrow\right\rangle_3}+\sqrt{\frac13}{\textstyle\left|1-1\right\rangle_{12}\left|\uparrow\right\rangle_3},
$$
$$
{\textstyle\left|\frac12\frac12(1)\frac12:\frac32-\frac32\right\rangle}=
{\textstyle\left|1-1\right\rangle_{12}\left|\downarrow\right\rangle_3}.
$$
Best Answer
There is no easy way to generalize this. In part this is because the basis states are not unique. In the case of 3 spin-$1/2$ particles, there are two sets of $S=1/2$ states and any linear combination of $\vert 1/2,1/2\rangle_1$ and $\vert 1/2,1/2\rangle_2$ is also a legitimate basis state.
The problem of multiplicities gets worse as you increase the number of spin-1/2 constituents. The number of times the final spin $S$ occurs in an $n$ spin-1/2 system is given in terms of dimensions of Young diagrams with $n$ boxes, and closely tied to the permutation of $n$ objects. This is well understood but constructing states is usually done algorithmically using a computer.
The "more" canonical approach is precisely to use the permutation group. In particular, using (for instance) class operators it's possible to construct different sets, but doing this manually is a serious pain. There's also an approach based on Young projection operators, but it also becomes rapidly impractical to do this "by hand".
You can check out
or
or even the venerable