The Hamiltonian is defined as \begin{equation} H = E(\mathbf{s}) \equiv - J\sum_i s_i s_{i+1} - h\sum_i s_i . \end{equation} where the state vector is $\mathbf{s}=(s_0,s_1,\ldots,s_{N-1})$ with $s_i = \pm 1$ for $i = 0\ldots N-1$ belonging to the least residue system modulo $N$, i.e., $\mathbb{Z}_N$, and corresponding to the periodic boundary condition.
The partition function is defined as \begin{equation} Z \equiv \sum_{\mathbf{s}} e^{-\beta E(\mathbf{s})}, \end{equation} which can be rewritten into a matrix form: \begin{align} Z & = \prod_i \left( \sum_{s_i=\pm 1} \right) e^{\sum_i \beta s_i(J s_{i+1} + h)} \notag \\ & = \prod_i \left( \sum_{s_i=\pm 1} e^{\beta s_i(J s_{i+1} + h)} \right) \notag \\ & = \operatorname{tr} \left(T^N\right) \end{align} where \begin{equation} T = \begin{pmatrix} e^{\beta(J-h)} & e^{\beta(-J-h)} \\ e^{\beta(-J+h)} & e^{\beta(J+h)} \end{pmatrix} \end{equation} is the transfer matrix.
The characteristic polynomial is \begin{equation} \lambda^2 - \left(e^{\beta (J+h)}+e^{\beta (J-h)}\right) \lambda + e^{2\beta J} - e^{-2\beta J}. \end{equation} Since the discriminant \begin{align} \mathrm{Disc}_\lambda & = \left(e^{\beta (J+h)}+e^{\beta (J-h)}\right)^2 - 4 \left(e^{2\beta J} - e^{-2\beta J}\right) \notag \\ & = e^{2\beta(J+h)} + e^{2\beta(J-h)} + 2 e^{2\beta J} - 4 \left(e^{2\beta J} - e^{-2\beta J}\right) \notag \\ & = \left(e^{\beta (J+h)}-e^{\beta (J-h)}\right)^2 + 4 e^{-2\beta J} \notag \\ & > 0, \end{align} the characteristic polynomial has two real roots, that is, the transfer matrix $T$ has two real eigenvalues.
The two eigenvalues are \begin{align} \lambda_\pm &= \frac{1}{2} \left[ \left( e^{\beta(J+h)}+e^{\beta(J−h)} \right) \pm \sqrt{ \left(e^{\beta (J+h)}-e^{\beta (J-h)}\right)^2 + 4 e^{-2\beta J} } \right] \notag \\ &= e^{\beta J}\left[\cosh(\beta h)\pm\sqrt{\sinh^2(\beta h)+e^{-4\beta J}}\right] \end{align} and the partition function is given by \begin{equation} Z = \lambda_+^N+\lambda_-^N . \end{equation}
The eigenvalues become \begin{equation} \lambda_\pm = e^{\beta}\pm e^{-\beta}. \end{equation} The partition function is \begin{equation} Z = \left(e^\beta+e^{-\beta}\right)^N+\left(e^\beta-e^{-\beta}\right)^N . \end{equation}
Generally, the average energy $U = \langle E\rangle$ can be calculated from the $\beta$-derivative of $-\log Z$: \begin{align} U &= \langle E\rangle \notag \\ &= Z^{-1} \sum_{\mathbf{s}} E e^{-\beta E} \notag \\ &= - Z^{-1} \frac{\partial}{\partial\beta} Z \notag \\ &= - \frac{\partial}{\partial\beta}\log Z \\ &= \frac{N\left[\left(e^\beta+e^{-\beta}\right)^{N-1}\left(e^\beta-e^{-\beta}\right)+\left(e^\beta+e^{-\beta}\right)\left(e^\beta-e^{-\beta}\right)^{N-1}\right]}{\left(e^\beta+e^{-\beta}\right)^N+\left(e^\beta-e^{-\beta}\right)^N} \notag \\ &= N\frac{r+r^{N-1}}{1+r^N} \end{align} where $r = \left(e^\beta-e^{-\beta}\right)/\left(e^\beta+e^{-\beta}\right)<1$ is the ratio between the two eigenvalues.