Often in condensed matter physics literature, one encounters a Hamiltonian that goes something like : $$ H = \sum_{i=1}^{n} J_{i}\ S_{i}^{z} S_{i+1}^{z}, $$ where $J_{i}$ are the coupling constants, $S_{i}^{z}$ representing the spin operator $S^{z}$ acting at $i^{\text{th}}$ site and $n$ being the total number of sites (with perhaps a ring-like identification). Now, $S^{z}$ is a $2 \times 2$ matrix with the spin at any particular site being a $2 \times 1$ matrix or column vector.
If one knew the wavefunction like $|\psi \rangle = \dots \otimes | i \rangle \otimes | i+1 \rangle \otimes \cdots $, then it would be easy to simply apply the Hamiltonian and do whatever one wants. But suppose that I don't have that $| \psi \rangle$ and need to directly diagonalize $H$ to calculate its eigenvalues & eigenvectors (for e.g., while simulating such a model on a lattice on computer), how would I go about doing that?
How would one make sure that $S^{z}_{i}$ is different from $S^{z}_{i+1}$? The only way to do that is to operate them on some presumed (random) spin existing at their sites. But if one does that then the Hamiltonian is a summation of product of two $2 \times 1$ matrices ($[2 \times 2]\ *\ [2 \times 1] \to [2 \times 1]$ matrix) which doesn't make sense!
Could someone elaborate where I'm going wrong? Basically, I want to diagonalize a Hamiltonian of the form above on a computer (using, say, Mathematica - which has inbuilt functions to give the eigenvalues & eigenvectors) without knowing the eigenfunctions.
Edit : Upon further thought, by saying that the spin at any particular site is a $2 \times 1$ matrix, am I implicitly assuming that it is a non-tensor state? Even if that's true, how does it help give the correct dimension of the Hamiltonian to diagonalize it?