Dimension of Hamiltonian & Diagonalizability

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?

Dimension of the Hilbert space

The Hilbert space of a spin chain of length $$n$$ is given by the tensor product of $$n$$ copies of the Hilbert space of a single spin. Assuming we are talking about spin-$$\frac{1}{2}$$, the dimension of the space is $$2^n$$.

In this context, an "operator acting on a single spin" means a $$2^n$$-dimensional operator obtained as the tensor product of a 2d operator with $$n-1$$ copies of the 2d identity, so that the action on the other spins is trivial.

For example, given the ($$2\times2$$) Pauli matrix $$\sigma^z$$,

$$S_i^z = \underbrace{\mathbb{I}_2 \otimes \cdots\otimes \mathbb{I}_2}_{i-1}\otimes \sigma^z\otimes\underbrace{\mathbb{I}_2\otimes \cdots\otimes \mathbb{I}_2}_{n-i} .$$

Similarly, a "two-body" operator $$S_i^z S_{i+1}^z$$ can be explicitly written

$$S_i^z S_{i+1}^z = \underbrace{\mathbb{I}_2 \otimes \cdots\otimes \mathbb{I}_2}_{i-1}\otimes \sigma^z\otimes \sigma^z\otimes\underbrace{\mathbb{I}_2\otimes \cdots\otimes \mathbb{I}_2}_{n-i-1},$$ and so on.

Numerical diagonalisation

Diagonalising such a Hamiltonian numerically is a hard problem. In principle, with the expressions above you can construct explicitly the matrix in your favourite programming language (e.g. in Mathematica using KroneckerProduct and PauliMatrix). However, as noted above, the dimension grows exponentially with the number of spins, meaning you will not be able even get to $$n\approx 20$$.

Therefore one needs to find smarter ways to do things; for example a starting point would be to represent $$H$$ as a sparse matrix. This is an active field of research, and even with state of the art algorithms and big computers, people are able to solve the problem with at most a few tens spins.