Energy Spectrum of pair of spin-1/2 particles with general Hamiltonian I found this problem, and so far I am stumped. I was wondering if anyone wanted to solve it with me, or help me calculate eigenvectors, or just give insight on my questions.

Consider a system of two spin-1/2 particles interacting through the Hamitonian
  $$H = A(S_x^2 − S_y^2) + BS_z^2,$$ 
  where $A$ and $B$ are constants and $S_x$, $S_y$ and $S_z$ are the three components of the total spin of the system. Find the energy spectrum and the corresponding eigenvectors.

It's been awhile since I did QM, and I know you can't have simultaneous eigenvectors of more than one measurement of spin along an axis, but how does that relate to this problem where there Sx^2 and Sy^2 are in the hamiltonian? I've been trying to work out 4x4 matrices, but I can't find good eigenvectors, but even then I am confused how applying a pauli matrix (say for x or y) transforms the eigenvector to another eigenvector, which I think is another example of not being able to define simultaneous eigenvectors etc. How would one start this problem? I will break it down to S^2 and Sz terms, but there will be remaining terms, will they just not contribute to the energy spectrum? If someone much more knowledgable than me could solve this problem easily and post it, that would clarify so much to me. For now, I'll try to read Feynman's lecture on the subject.
 A: It's not clear to me at this point that there is a tricky way to do this problem, but here is a systematic way that would work for any hamiltonian written in terms of the components of the total spin operator.
As you have indicated, the Hamiltonian acts on a $4$-dimensional Hilbert space $\mathcal H$, so it can be written as a $4$-by-$4$ matrix in a given basis.  The eigenvalues of the matrix don't depend on the basis you choose, so provided you can determine the matrix representation of $H$ in some basis, you can just find the eigenvalues of this matrix, and you'll be done.
Is there a convenient basis to choose?  Well, I'd say that the basis $|s, m\rangle$ consisting of simultaneous eigenvectors of $\mathbf S^2$ and $S_z$, the total spin squared and total $z$-component of spin operators, is best because we'll easily be able to compute the matrix elements of the Hamiltonian in this basis as you'll see in a moment;
$$
  \mathbf S^2|s,m\rangle = \hbar^2s(s+1)|s,m\rangle, \qquad S_z |s,m\rangle = \hbar m|s,m\rangle
$$
Since the system consists of two spin-$1/2$ particles, this basis contains three spin-$1$ states and one spin-$0$ state; $\{|1,1\rangle,|1,0\rangle, |1,-1\rangle, |0,0\rangle\}$ .  These are the so-called triplet and singlet states that arise from appropriately "adding" the two spins quantum mechanically.
Now, of course the $S_z^2$ term is particularly easy to deal with in this basis, in fact, the matrix representation of $S_z^2$ in this basis is diagonal.  The trick is to therefore compute the matrix representation of $S_x^2 - S_y^2$ in this basis.  This can be done systematically using raising and lowering operators.  Recall that
$$
  S_+ =S_x + iS_y, \qquad S_- = S_x - iS_y
$$
so that
$$
  S_x = \frac{1}{2}(S_+ + S_-), \qquad S_y = \frac{1}{2i}(S_+ - S_-) 
$$
and therefore for example
$$
  S_x^2 = \frac{1}{4}(S_+^2 + 2S_+S_- + S_-^2)
$$
and similarly for $S_y^2$.  Now you just need to use these expressions to determine the matrix representation of the Hamiltonian in this basis and find its eigenvalues.
I strongly suspect that there is a slicker way of doing this problem, but it's just not coming to me at the moment.
A: 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}$.
A: If one is able to have an explicit 4 * 4 matrix for the hamiltonian, the work to find eigenvalues and eigenvectors is a basic math problem. So I will suppose you know that.
I will answer to how get this hamiltonian matrix.
1) a basis, for the 2-spin state is $|0>\otimes|0>, ~|0>\otimes|1>, ~|1>\otimes|0>, ~|1>\otimes|1>$
2) The meaning of $S_i$ is in fact $S_i = ((s_1)_i \otimes \mathbb{Id_2} + \mathbb{Id_1} \otimes (s_2)_i)$, 
where $(s_1)_i$ is a  $2 * 2$ matrix (spin operator) acting on the first particle, and $(s_2)_i$ is a $2 * 2$ matrix (spin operator) acting on the second particle, and where $\mathbb{Id_1}$ and $\mathbb{Id_2}$ are $2 * 2$ matrix acting respectively on particules #1 and #2.
Note, then you have $(s_1)_i = (s_2)_i = \frac{1}{2} \sigma_i$, where  $\sigma_i$ is the Pauli matrice.
Note, then, that you have  $((s_1)_i)^2 = ((s_2)_i)^2 = \frac{1}{4} \mathbb{Id}$
3) Now, you can simplify the hamiltonian : $H = A(S_x^2 − S_y^2) + BS_z^2,$
For instance $S_x^2 = ((s_1)_x \otimes \mathbb{Id_2} + \mathbb{Id_1} \otimes (s_2)_x)^2 = \frac{1}{2} \mathbb{Id_1} \otimes \mathbb{Id_2} + \frac{1}{2}\sigma_x \otimes \sigma_x$.
So, you will have : 
$H =\frac{B}{2}  \mathbb{Id} + \frac{1}{2} (A(\sigma_x \otimes \sigma_x - \sigma_y \otimes \sigma_y) + B \sigma_z \otimes \sigma_z)$
Here $\mathbb{Id} = \mathbb{Id_1} \otimes \mathbb{Id_2}$ is the $4*4$  identity matrix.
4) To get the hamiltonian matrix, now just apply H to the basic vectors, for instance : 
$$H ~(|0>\otimes|0>) = \frac{B}{2} |0>\otimes|0> + \frac{1}{2}(A(|1>\otimes|1> - i^2 |1>\otimes|1>) 
+ B|0>\otimes|0>)$$
That is : $$H ~(|0>\otimes|0>) = B|0>\otimes|0> + A |1>\otimes|1>$$
By doing the same thing with the other basis vectors, you will get the full  $4*4$ hamiltonian matrix.
Note: if you understand the notation $|0>\otimes|0>$, you can use the notation $|0>|0>$ which is simpler.
