I've tried an alternative solution for finding the partition function of this model.
So is what I've done correct? If it isn't then please prove and explain why not. (I'm pretty sure I made a mistake along the way):$$ $$
Starting with the Energy:
$$
E = -J\sum_{i=1}^{N}S_{i}S_{i+1}-J'\sum_{i=1}^{N}S_{i}S_{i+2}-B\sum_{i=1}^{N}S_i
$$
I've used the same trick used for NN but instead of two terms for the energy contribution from B I now have 3:
$$
E = -J\sum_{i=1}^{N}S_{i}S_{i+1}-J'\sum_{i=1}^{N}S_{i}S_{i+2}-\frac{B}{3}\sum_{i=1}^{N}(S_i+S_{i+1}+S_{i+2})
$$
The point of doing this becomes more obvious when I write out the partition function:
$$
Z=\sum_{\{S_{i}\}}e^{-\beta H}
$$
$$
=\sum_{S_i=\pm1}…\sum_{S_N=\pm1}e^{\beta\sum_{i}[JS_{i}S_{i+1}+J'S_{i}S_{i+2}+\frac{B}{3}(S_i+S_{i+1}+S_{i+2})]}
$$
Having written the energy in a way that I think is analogous to Isings method allows me to write this as the product of two exponentials and define two transfer matrices as follows:
$$
Z=\sum_{S_i=\pm1}…\sum_{S_N=\pm1}e^{\beta\sum_{i}(JS_{i}S_{i+1}+\frac{B}{3}(\frac{S_i}{2}+S_{i+1}))}e^{\beta\sum_{i}(J'S_{i}S_{i+2}+\frac{B}{3}(\frac{S_i}{2}+S_{i+2}))}
$$
And introducing $T_{_{NN}}= \left( \begin{array}{cc}
e^{\beta(J+\frac{B}{2})} & e^{-\beta(J+\frac{B}{6})} \\
e^{\beta(-J+\frac{B}{6})} & e^{\beta(J-\frac{B}{2})} \\
\end{array} \right) $ and $T_{_{NNN}}= \left( \begin{array}{cc}
e^{\beta(J'+\frac{B}{2})} & e^{-\beta(J'+\frac{B}{6})} \\
e^{\beta(-J'+\frac{B}{6})} & e^{\beta(J'-\frac{B}{2})} \\
\end{array} \right)$
Writting out the partition function in terms of these two transfer matrices and rearranging:
$$Z=\sum_{S_1=\pm 1}…\sum_{S_N=\pm 1}\langle S_1| T_{_{NN}}|S_2\rangle\langle S_1| T_{_{NNN}}|S_3\rangle\langle S_2| T_{_{NN}}| S_3\rangle\langle S_2| T_{_{NNN}}|S_4\rangle…\langle S_N| T_{_{NN}}|S_1\rangle\langle S_N| T_{_{NNN}}|S_2\rangle$$
$$=\sum_{S_1=\pm 1}…\sum_{S_N=\pm 1}\langle S_2| T_{_{NN}}^{\dagger}|S_1\rangle\langle S_1| T_{_{NNN}}|S_3\rangle\langle S_3| T_{_{NN}}^{\dagger}| S_2\rangle\langle S_2| T_{_{NNN}}|S_4\rangle…\langle S_1| T_{_{NN}}^{\dagger}|S_N\rangle\langle S_N| T_{_{NNN}}|S_2\rangle$$
And using completeness :
$$=\sum_{S_2=\pm1}\langle S_2|( T_{_{NN}}^{\dagger}T_{_{NNN}})^N|S_2\rangle$$
$$=Tr( T_{_{NN}}^{\dagger}T_{_{NNN}})^N$$
And voila, I have something I can plug into matlab and calculate U,F,M,$\chi$, heat capacity etc. for arbitrary J,J',B T and N (which I've done). The only problem is that it doesn't give me the same results as the NN 1D ising for J'=0 and it doesn't give me exactly the same results as the 4X4 matrix shown here: https://arxiv.org/abs/cond-mat/9703187 (equation 9).
So I must have made some kind of mistake along the way, however, I would like to point out that the 4X4 matrix doesn't seem to reduce to the 1D ising model with NN interaction when J'=0 either.
(I should specify that my matlab numerical computations have so far used values of N ranging from 20 to 500, I've also tried calculating U/N and $\frac{\chi}{N}$using the largest eigenvalues from each transfer matrix for certain parameters and obtained the same results for T>0)
Best Answer
Note that in the sum
$$=\sum_{S_1=\pm 1}...\sum_{S_N=\pm 1}\langle S_2| T_{_{NN}}^{\dagger}|S_1\rangle\langle S_1| T_{_{NNN}}|S_3\rangle\langle S_3| T_{_{NN}}^{\dagger}| S_2\rangle\langle S_2| T_{_{NNN}}|S_4\rangle...\langle S_1| T_{_{NN}}^{\dagger}|S_N\rangle\langle S_N| T_{_{NNN}}|S_2\rangle$$
every pair $|S_i\rangle\langle S_i|$ occurs twice with some operator/matrix between them. So you cannot simply execute the sum over $S_i = \pm1$ for both occurrences separately.
As pointed out in the paper you linked, one solution via transfer matrix is to group two neighbouring spins to one 4-state spin. Then the problem can be reduced to NN interactions only.