As mentioned in the comments. Your matrix representation of the creation and annihilation operators is incorrect. This is easy to see since
\begin{align}
a ^\dagger _{ nm} & = \left\langle n \right| a ^\dagger \left| m \right\rangle
\\
& = \sqrt{ m + 1 } \delta _{ n , m + 1 }.
\end{align}
Thus we have,
\begin{equation}
\left( \begin{array}{cccc}
0 & 0 & 0 & ...\\
1 & 0 & 0 &...\\
0 & \sqrt{2} & 0 &... \\
0 &0 & \sqrt{3} & \ddots\\
\vdots &... & \ddots & \ddots
\end{array} \right)
\end{equation}
The kets in this space are simple:
\begin{equation}
\left| n \right\rangle = \left( \begin{array}{c}
0 \\
\vdots \\
0 \\
1 \\
0 \\
\vdots \\
0
\end{array} \right)
\end{equation}
where $1$ is at the $ n +1 $'th spot down the vector. Its easy to see that acting on this vector we get a result consistent with the relation for the raising operator above.
The vacuum state is then just given by
\begin{equation}
\left| 0 \right\rangle = \left( \begin{array}{c}
1 \\
0 \\
\vdots \\
0
\end{array} \right)
\end{equation}
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.
Best Answer
First of all, the equation $$ \begin{equation}A\otimes B=A\otimes \mathbb{1}+\mathbb{1}\otimes B,\end{equation} $$ is a claim about an identity, and this claim is incorrect. Note that for $1\times 1$ matrices, the matrices are numbers and the equation above reduces to $$ a\cdot b = a\cdot 1 + 1 \cdot b$$ which is clearly wrong because the addition (the right hand side is just $a+b$) is something else than multiplication!
If we have a Hilbert space $H = H_A \otimes H_B$ and there is a change of the basis in the Hilbert space $H_A$ separately and in $H_B$ separately, the change of the basis in the Hilbert space $H = H_A \otimes H_B$ is the tensor product of the two transformation matrices, $A\otimes B$, in your notation (also expressed in the block form using those $a_{11}B$ etc.), and not $A\otimes 1 + 1\otimes B$.
This is easily calculated. The basis vectors transform as $$ a_j = a'_i \cdot A_{ij}, \quad b_j = b'_i\cdot B_{ij} $$ (or switch the order of the indices or which side has the prime or invert the matrix or transpose it etc. – at least some of these choices are conventions; you must be careful to use consistent conventions all the time) and the basis vectors of $H_A\otimes H_B$ are $$ a_j\otimes b_\ell = a'_i \cdot A_{ij}\otimes b'_k\cdot B_{k\ell} = (a'_i\otimes b'_k) \cdot A_{ij}B_{k\ell} $$ The rule here is that algebraically, the tensor multiplication is treated just like a normal multiplication (with no identification of indices or summing over indices). The transformation matrix in $H_A\otimes H_B$ is therefore $A_{ij}B_{k\ell}$ where $i,k$ are parts of the "first" index generalizing $i$ or $k$, and $j,\ell$ are analogously the second index.
The expression $$A\otimes 1 + 1\otimes B$$ is also important and it appears at many places, but we must discuss the "derivatives", the changes under the infinitesimal transformations. It's because the expression above is analogous to the Leibniz rule for the derivative of the product $$ (ab)' = a'b + ab' $$ In quantum mechanics, derivatives may also be represented as commutators. For example, if $\vec J$ is the operator of the total angular momentum, then $$ [J_z,A] $$ is the operator measuring the $z$-component of the angular momentum carried by the operator $A$. Similarly for $B$. Because we may decompose the commutator with a product to the sum of two commutator-like terms, $$ [J_z,AB] = J_z AB - AB J_z = J_z AB - A J_z B + A J_z B - AB J_z = [J_z,A]B+A[J_z,B]$$ one may see that the result of acting with $J_z$ on the whole $AB$ is a sum of two terms – the angular momentum of $A$ and the angular momentum of $B$. Those give the two terms you mentioned because the angular momentum of $A$ only acts as $J_z^A$ on the space $H_A$, but it acts as the identity operator (one) on the space $H_B$, and vice versa, and that is why you have the tensor products with the identity matrix.
So the expressions $A\otimes 1+1\otimes B$ would appear but not if you considered a "finite" change of the bases, but if you only considered changes of bases that are "infinitesimally close" to the identity change (no change at all), and if you would take a difference between these two expressions (a derivative with respect to some parameters labeling the change of the bases, taken near the trivial transformation).
In terms of groups and algebras, the straightforward $A\otimes B$ gives us the matrix of a finite transformation in the Lie group as it acts on the tensor-product representation $H_A\otimes H_B$. On the other hand, $A\otimes 1+1\otimes B$ is the form of the Lie algebra generator with respect to the same space. Elements of Lie algebras are matrices/operators $(g-1)/\epsilon$ where $g$ is a group element infinitesimally close to the identity $1$ of the Lie group.