Hamiltonian Matrix for XXZ Model Given the XXZ model Hamiltonian,
$H = -\frac{1}{2}\sum^{N}_{i}(\sigma_{i}^{x}\sigma_{i+1}^{x}+\sigma_{i}^{y}\sigma_{i+1}^{y}+\Delta\sigma_{i}^{z}\sigma_{i+1}^{z})$
The two-site Hamiltonian reads
$H = \begin{bmatrix}-\Delta&0&0&0\\0&\Delta&-2&0\\0&-2&\Delta&0\\0&0&0&-\Delta\end{bmatrix}$
How do I get the Hamiltonian matrix from the Hamiltonian?  
 A: In order to fix notation, for a spin chain with $N$ places/sites, the two-sites hamiltonian is 
$$
H_{i,i+1}=-\frac{1}{2}\left(J_x\sigma^x_{i}\otimes\sigma^x_{i+1}+J_y\sigma^y_{i}\otimes\sigma^y_{i+1}+J_z\sigma^z_{i}\otimes\sigma^z_{i+1}\right),
$$
 for $J_x,J_y,J_z$ are three scalars, not necesarily different, so the hamiltonian for the complete chain $H=\sum_{i=1}^N H_{i,i+1}$, where $i+N\equiv i$ in the subindices (imposing in this way periodic boundary conditions in the chain). 
In each site of the chain a state is a vector of $Span(|\uparrow\rangle,|\downarrow\rangle)\simeq\mathbb{C}^2$ (so for the full hamiltonian the Hilbert space would be the tensor product of the hilbert spaces in each site $\bigotimes_{i=1}^N\mathbb{C}^2$).
Focusing in the two-site hamiltonian, recall how it acts on a space. A state of two spins will be $|\psi_1\rangle\otimes|\psi_2\rangle$, where $|\psi_i\rangle$ will be a certain linear combination of $|\uparrow\rangle $ and $|\downarrow\rangle$, so for an operator $A_1B_2$
$$A_1B_2|\psi_1\rangle\otimes|\psi_2\rangle=(A_1|\psi_1\rangle)\otimes(B_2|\psi_2\rangle).$$
To symplify notation, I will get rid of the symbol $\otimes$, keeping the subindices as an implicit way to show how the operators act on states. To check if the hamiltonian is (skew-)symmetric, you need to show that $H_{1,2}=H_{1,2}^\dagger$, or equivalently, $\forall|\psi_1\rangle|\psi_2\rangle,|\phi_1\rangle|\phi_2\rangle,$
$$\langle\phi_1|\langle\phi_2| H_{1,2}|\psi_1\rangle|\psi_2\rangle=(\langle\psi_1|\langle\psi_2|H_{1,2}|\phi_1\rangle|\phi_2\rangle)^\dagger.$$ 
Here it is shown for one term, but the generalization is straightforward:
$$\langle\phi_1|\langle\phi_2| \sigma^x_1\sigma^x_2|\psi_1\rangle|\psi_2\rangle=\\
=\langle\phi_1| \sigma^x_1|\psi_1\rangle\langle\phi_2|\sigma^x_2|\psi_2\rangle=\\
=(\langle\psi_1| {\sigma^x_1}^\dagger|\phi_1\rangle)^\dagger(\langle\psi_2|{\sigma^x_2}^\dagger|\phi_2\rangle)^\dagger=\\
=(\langle\psi_1| {\sigma^x_1}|\phi_1\rangle)^\dagger(\langle\psi_2|{\sigma^x_2}|\phi_2\rangle)^\dagger=\\
=(\langle\psi_1| {\sigma^x_1}|\phi_1\rangle\langle\psi_2|{\sigma^x_2}|\phi_2\rangle)^\dagger=\\
=(\langle\psi_1|\langle\psi_2| {\sigma^x_1}{\sigma^x_2}|\phi_1\rangle|\phi_2\rangle)^\dagger,
$$
where in the third to fourth step I have used $\sigma=\sigma^\dagger$ for Pauli matrices. The term is (skew-)symmetric, and the hamiltonian is sum of similar terms, so it must be (skew-)symmetric as well. 
If you choose a specific basis of states for the Hilbert space, you can get the matricial form of the hamiltonian. I use the same term as before to work and as a bases of states $\{|\uparrow\rangle|\uparrow\rangle,|\uparrow\rangle|\downarrow\rangle,|\downarrow\rangle|\uparrow\rangle, |\downarrow\rangle|\downarrow\rangle\}$. Recall that:
$$\sigma^x|\uparrow\rangle=|\downarrow\rangle,\quad \sigma^x|\downarrow\rangle=|\uparrow\rangle,$$
so the matrix representation  of $\sigma^x$ is given by $\left(\matrix{0 & 1 \\ 1 & 0}\right)$, and then, the matrix representation of the tensor product:
$$
\sigma_1^x\sigma_2^x=\left(\matrix{0 & 1 \\ 1 & 0}\right)\otimes\left(\matrix{0 & 1 \\ 1 & 0}\right)=\left(\matrix{0\cdot\left(\matrix{0 & 1 \\ 1 & 0}\right) & 1\cdot\left(\matrix{0 & 1 \\ 1 & 0}\right) \\ 1\cdot\left(\matrix{0 & 1 \\ 1 & 0}\right) & 0\cdot\left(\matrix{0 & 1 \\ 1 & 0}\right)}\right)=\left(\matrix{0 & 0 & 0 & 1 \\0 & 0 & 1 & 0 \\0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0}\right)
$$
For the a general $J_x,J_y,J_z$
$$
H_{i,i+1}=\left(\matrix{J_z & 0 & 0 & J_x-J_y \\0 & -J_z & J_x+J_y & 0 \\0 & J_x+J_y & -J_z & 0 \\ J_x-J_y & 0 & 0 & J_z}\right)
$$
