Here I basically do what joshphysics has already mentioned, just in a little more detail, and in a bit more intuitive basis (which makes effectively no difference). So definitely not a slick way.
I use the basis $|m_1\rangle\otimes|m_2\rangle$, where $m_i=\pm1/2$. Keeping only the sign we denote the basis as $\{|++\rangle,|+-\rangle,|-+\rangle,|--\rangle\}$.
Now for $S_z=S_{1z}+S_{2z}$ you can check that
$$
S_z|++\rangle=\left(\frac{\hbar}{2}+\frac{\hbar}{2}\right)|++\rangle
$$
and similarly
$$
S_z|--\rangle=\left(-\frac{\hbar}{2}-\frac{\hbar}{2}\right)|--\rangle
$$
so that you and up with the matrix
$$
S_z\rightarrow\hbar\left[\begin{array}{cccc}1&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&-1\end{array}\right]
$$
and
$$
S_z^2\rightarrow\hbar^2\left[\begin{array}{cccc}1&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&1\end{array}\right]
$$
For $S_x$ and $S_y$ you use raising and lowering operators as josh suggested and obtain the matrix elements by acting on each state in the basis and seeing what has a non-zero dot product with the result, e.g.
$$
\begin{align}
S_y|+-\rangle=\frac{1}{2i}\left(S_+-S_-\right)|+-\rangle&=\\\frac{1}{2i}\left(S_{1+}+S_{2+}-S_{1-}-S_{2-}\right)|+-\rangle&=\frac{1}{2i}\left(\hbar|++\rangle-\hbar|--\rangle\right)
\end{align}
$$
thus $\langle++|S_y|+-\rangle=\frac{\hbar}{2i}$ and $\langle--|S_y|+-\rangle=-\frac{\hbar}{2i}$. Once you do this you can get
$$
S_x\rightarrow\hbar\left[\begin{array}{cccc}0&1&1&0\\1&0&0&1\\1&0&0&1\\0&1&1&0\end{array}\right]
$$
and
$$
S_y\rightarrow\hbar\left[\begin{array}{cccc}0&1&1&0\\-1&0&0&1\\-1&0&0&1\\0&-1&-1&0\end{array}\right]
$$
Finally the Hamiltonian looks like
$$
H=\left[\begin{array}{cccc}\hbar^2B&0&0&\hbar^2A\\0&0&0&0\\0&0&0&0\\\hbar^2A&0&0&\hbar^2B\end{array}\right]
$$
Therefore $|+-\rangle$ and $|-+\rangle$ are eigenvectors with eigenvalue $0$, and you only need diagonalize $2\times2$ matrix to find that the eigenvalue $\hbar^2\left(B+A\right)$ corresponds to the eigenvector $\left(|++\rangle+|--\rangle\right)/\sqrt{2}$, and $\hbar^2\left(B-A\right)$ corresponds to the eigenvector $\left(|++\rangle-|--\rangle\right)/\sqrt{2}$.
The two particles $m_s$ and $m_I$ live in different vector spaces, so you are actually not picking the same basis vectors (because the basis vectors of the different particles belong to two separate vector spaces).
Secondly, the tensor product between the basis vectors of the two different vector spaces will form the basis vectors of a new $3 \times 3 = 9$ dimensional vector space. For instance:
\begin{equation}
|m_s=1\rangle \otimes |m_I=0\rangle = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}
\end{equation}
So the equation you have written down is correct.
I would strongly recommend these short set of notes to get a better feeling of the full machinery.
Edit in response to comment:
We need to use the tensor product to add angular momenta. For simplicity, consider a system with two particles, and the spin of both particles is $1/2$. Then, these particles can be combined in the following way:
\begin{equation}
\begin{array}{cccc}
\uparrow \uparrow \; ,& \uparrow \downarrow \; ,& \downarrow \uparrow \; ,& \downarrow \downarrow
\end{array}
\end{equation}
This means that the two-particle Hilbert space (i.e. the Hilbert space corresponding to the system) is spanned by four basis vectors:
\begin{equation}
|s_1,m_1; s_2,m_2 \rangle \equiv | s_1,m_1 \rangle \otimes |s_2,m_2\rangle
\end{equation}
Subsequently, as I have sketched in this answer, we decompose this tensor product to determine what the possible eigenvalues for the magnitude and $z$-component of the total system (i.e. the system formed by the two spin-$1/2$ particles) can be.
Best Answer
Hints:
Notice that if we define $\mathbf S_{123} = \mathbf S_1 + \mathbf S_2 + \mathbf S_3$ and $\mathbf S_{12} = \mathbf S_1 + \mathbf S_2$, then we have \begin{align} \mathbf S_{123}^2 = \mathbf S_{12}^2 + \mathbf S_3^2 + 2\mathbf S_{12}\cdot\mathbf S_3 \end{align}
Notice that your hamiltonian can be written as follows: \begin{align} H = \frac{A}{2}(\mathbf S_{12}^2- \mathbf S_1^2 - \mathbf S_2^2) + B(\mathbf S_{12}\cdot\mathbf S_3) \end{align}
Combine the last two hints.
Recall that the representation theory of the angular momentum algebra (aka addition of angular momentum) tells us that a tensor product of three spin-$1/2$ reprsentations splits into a direction sum as follows: \begin{align} \tfrac{1}{2}\otimes\tfrac{1}{2}\otimes\tfrac{1}{2} &= (1 \oplus 0)\otimes\tfrac{1}{2} \\ &= (1\otimes \tfrac{1}{2}) \oplus (0\otimes \tfrac{1}{2}) \\ &= \underbrace{(\tfrac{3}{2}\oplus\tfrac{1}{2})}_{s_{12}=1, s_3 = \frac{1}{2}}\oplus \underbrace{\tfrac{1}{2}}_{s_{12}=0, s_3=\frac{1}{2}} \end{align}
(Addendum) To clarify some discussion in the comments, the last line in number 4 means that the Hilbert space admits an orthonormal basis of states which I'll label $|s_{123}, s_{12}\rangle$ for which \begin{align} \mathbf S_{123}^2|s_{123}, m_{123},s_{12}\rangle &= \hbar^2 s_{123}(s_{123}+1)|s_{123}, m_{123},s_{12}\rangle \\ S_{123}^z |s_{123}, m_{123},s_{12}\rangle &= \hbar m_{123}|s_{123}, m_{123},s_{12}\rangle\\ \mathbf S_{12}^2|s_{123}, m_{123},s_{12}\rangle &= \hbar^2 s_{12}(s_{12}+1)|s_{123}, m_{123},s_{12}\rangle \\ \mathbf S_{1}^2|s_{123}, m_{123},s_{12}\rangle &= \hbar^2 \tfrac{1}{2}(\tfrac{1}{2}+1)|s_{123}, m_{123},s_{12}\rangle \\ \mathbf S_{2}^2|s_{123}, m_{123},s_{12}\rangle &= \hbar^2 \tfrac{1}{2}(\tfrac{1}{2}+1)|s_{123}, m_{123},s_{12}\rangle \\ \mathbf S_{3}^2|s_{123}, m_{123},s_{12}\rangle &= \hbar^2 \tfrac{1}{2}(\tfrac{1}{2}+1)|s_{123}, m_{123},s_{12}\rangle \end{align} and explicitly, the basis is as follows: \begin{align} \left.\begin{array}{l} |\tfrac{3}{2}, \tfrac{3}{2},1\rangle \\ |\tfrac{3}{2}, \tfrac{1}{2},1\rangle \\ |\tfrac{3}{2}, -\tfrac{1}{2},1\rangle \\ |\tfrac{3}{2}, -\tfrac{3}{2},1\rangle \end{array}\right\} s_{123} = \tfrac{3}{2}, s_{12} = 1, \dim = 4\\ \left.\begin{array}{l} |\tfrac{1}{2}, \tfrac{1}{2},1\rangle \\ |\tfrac{1}{2}, -\tfrac{1}{2},1\rangle \end{array}\right\} s_{123} = \tfrac{1}{2}, s_{12} = 1, \dim = 2\\ \left.\begin{array}{l} |\tfrac{1}{2}, \tfrac{1}{2},0\rangle \\ |\tfrac{1}{2}, -\tfrac{1}{2},0\rangle \end{array}\right\} s_{123} = \tfrac{1}{2}, s_{12} = 0, \dim = 2 \end{align} This is a basis for the entire three spin-$1/2$ particle Hilbert space; it is eight-dimensional. There are therefore a total of 8 orthonormal eigenvectors for the hamiltonian.