I am trying to derive the spin matrices $S_x$ and $S_y$ in the z-basis for spin $s>1/2$ using the direct product (Kronecker product) method. For simplicity, let's focus on the case $s=1$.
I understand how to derive $S_z$. As shown in https://physics.stackexchange.com/a/342156/50484,
if for example
\begin{equation}
S_{1z}=\tfrac{1}{2}
\begin{bmatrix}
1 & 0\\
0 &\!\!\! -\!1
\end{bmatrix}\;,\;
S_{2z}=\tfrac{1}{2}
\begin{bmatrix}
1 & 0\\
0 &\!\!\! -\!1
\end{bmatrix}
\tag{08}
\end{equation}
then
\begin{equation}
S_{z-tot}=\left(S_{1z} \otimes I_2\right)+ \left(I_1 \otimes S_{2z}\right)=
\begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0\\
0 & 0 & 0 &\!\!\! -\!1
\end{bmatrix}
\tag{09}
\end{equation}
The matrix in (09) is already diagonal with eigenvalues 1,0,0,-1. Rearranging rows and columns we have
\begin{equation}
S'_{z-tot}=
\begin{bmatrix}
\begin{array}{c|cccc}
0 & 0 & 0 & 0\\
\hline
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0\\
0 & 0 & 0 & \!\!\!\!-\!1
\end{array}
\end{bmatrix}
=
\begin{bmatrix}
\begin{array}{c|c}
S_{z}^{(j=0)} & 0_{1\times 3}\\
\hline
0_{3\times1} & S_{z}^{(j=1)}
\end{array}
\end{bmatrix}
\tag{10}
\end{equation}
and so
$S_{z,s=1}=
\begin{bmatrix}
1 & 0 & 0\\
0 & 0 & 0 \\
0 & 0 & -1\\
\end{bmatrix}$.
It is precisely the steps starting from (09) and getting to $S_{z,s=1}$ that I am wondering about, but for $S_{x,s=1}$ and $S_{y,s=1}$. Specifically, I can get
\begin{equation}
S_{x-tot}=\left(S_{1x} \otimes I_2\right)+ \left(I_1 \otimes S_{2x}\right)
=\tfrac{1}{2}
\begin{bmatrix}
0 & 1 & 1 & 0\\
1 & 0 & 0 & 1 \\
1 & 0 & 0 & 1\\
0 & 1 & 1 & 0
\end{bmatrix}
\tag{20}
\end{equation}
\begin{equation}
S_{y-tot}=\left(S_{1y} \otimes I_2\right)+ \left(I_1 \otimes S_{2y}\right)=\tfrac{1}{2}
\begin{bmatrix}
0 & \!\!\! -\!i & \!\!\! -\!i & 0\\
i & 0 & 0 & \!\!\! -\!i \\
i & 0 & 0 & \!\!\! -\!i \\
0 & i & i & 0
\end{bmatrix}
\tag{21}
\end{equation}
but from these I do not know how to get
$S_{x,s=1}= \frac{1}{\sqrt 2}
\begin{bmatrix}
0 & 1 & 0\\
1 & 0 & 1 \\
0 & 1 & 0\\
\end{bmatrix}$
and
$S_{y,s=1}= \frac{1}{\sqrt 2}
\begin{bmatrix}
0 & -i & 0\\
i & 0 & -i \\
0 & i & 0\\
\end{bmatrix}$.
It seems $S_{z-tot}$ is just abnormally easy to work with because it is a diagonal matrix.
Best Answer
The problem with your approach is there are two $m=0$ states and the $S=1$ states are linear combination of those. In other words, you have hastily identified the $(1,1)$ as entirely in the $S=0$, and the $(3,3)$ as entirely in the $S=1$ subspace. In fact, it's a linear combination of those basis states which are in the relevant subspace.
Making the change of basis will solve your problem. You can make this change of basis "by hand" by searching for the linear combination of your first and third basis states that is killed by $S_+$, identifying this linear combination with the $S=0$ states.
Another method is to take the $m=1$ state (it is unique), and use $S_-=S_x-iS_y$ to lower to a single state, which must be the $S=1,m=0$ state.
A more elegant method uses Clebsch-Gordan coefficients, and another more computationally intensive method is to diagonalize $L^2=L_x^2+L_y^2+L_z^2$ to transform to the correct $S$ subspace.
Edit: For your single-particle states you are using the basis $\{\vert +\rangle,\vert -\rangle\}$ since the diagonal form of $S_z$ produces the eigenvalues $\pm \frac{1}{2}$ on the basis vectors $\vert+\rangle\mapsto (1,0)^\top$ and $\vert -\rangle \mapsto (0,1)^\top$.
A basis for your tensor product is thus $\{\vert +\rangle\vert +\rangle,\vert +\rangle\vert -\rangle,\vert -\rangle\vert +\rangle,\vert -\rangle\vert -\rangle$ with eigenvalues of your $S_{z-tot}$ given by $1,0,0,-1$ respectively.
The states $\vert +\rangle\vert +\rangle$ and $\vert -\rangle\vert -\rangle$ must live in the $S=1$ and must be $\vert 1,1\rangle$ and $\vert 1,-1\rangle$ by inspection. On the other hand the states $\vert +\rangle\vert -\rangle$ and $\vert +\rangle\vert -\rangle$ do not live in an irreducible subspace. This is immediately seen by lowering from $\vert +\rangle\vert +\rangle$.
Indeed, the $S=1,m=0$ state is the combination $$ \vert 1,0\rangle = C^{1,0}_{1/2,1/2;1/2,-1/2}\vert +\rangle \vert -\rangle +C^{1,0}_{1/2,-1/2;1/2,1/2}\vert -\rangle \vert +\rangle $$ where $C^{1,0}_{1/2,1/2;1/2,-1/2}$ is a Clebsch-Gordan coefficient. Likwise $$ \vert 0,0\rangle = C^{0,0}_{1/2,1/2;1/2,-1/2}\vert +\rangle \vert -\rangle +C^{0,0}_{1/2,-1/2;1/2,1/2}\vert -\rangle \vert +\rangle\, . $$ The general change of basis transformation $$ \vert S,m_s\rangle =\sum_{m_1m_2} C^{Sm}_{1/2,m_1;1/2,m_2}\vert m_1\rangle \vert m_2\rangle $$ defines a matrix $T$ with elements $C^{Sm}_{m_1m_2}$ that will bring your $4\times 4$ matrices to the correct block diagonal form.