Multiple summation for 1D Ising model

Question

I am trying to perform multiple summation in a Matlab code, but I don't know how to write a code to perform multiple for's . $$\sum_{S_1=\pm1}...\sum_{S_N=\pm1}(\Pi_{i=1}^{N-1}e^{S_iS_{i+1}})=\sum_{S_1=\pm1}...\sum_{S_N=\pm1}(e^{S_1S_2...S_{N-1}S_N})$$ As you can see if I try to implement this manually I should write N for's and 1 more for the big pi, so how can I reduce the N for's so I only need to use 3-4, one indicating the number of for that I want to repeat, another for the summation and the last one for the multiplication. If I do it manually for 10 spins, what I have to write is:

for S1 = [1, -1]
    for S2= [1, -1]
        a2 = S1*S2;
        for S3 = [1, -1]
            a3 = a2 + S2*S3;
            for S4 = [1, -1]
                a4 = a3+ S3*S4;
                for S5 = [1, -1]
                    a5 = a4 + S4*S5;
                    for S6 = [1, -1]
                        a6 = a5 + S5*S6;
                        for S7= [1, -1]
                            a7 = a6 + S6*S7;
                            for S8 = [1, -1]
                                a8 = a7 + S7*S8;
                                for S9 = [1, -1]
                                    a9 = a8 + S8*S9;
                                    for S10 = [1, -1]
                                        a10 = a9 + S9*S10;
                                        if S10 == 1
                                            zold10 = exp(a10);
                                        else
                                            z = exp(a10)+zold10;
                                        end
                                    end
                                    if S9 == 1
                                        zold9 = z;
                                    else 
                                        z = z + zold9;
                                    end
                                end
                                if S8 == 1
                                    zold8 = z;
                                else
                                    z = z + zold8;
                                end
                            end
                            if S7 == 1
                                zold7 = z; 
                            else
                                z = z+ zold7;
                            end
                        end
                        if S6 == 1
                            zold6 = z; 
                        else
                            z = z+ zold6;
                        end
                    end
                    if S5 == 1
                        zold5 = z;
                    else
                        z = z+zold5;
                    end
                end
                if S4 == 1
                    zold4 = z;
                else 
                    z = z + zold4;
                end
            end
            if S3 == 1
                zold3 = z;
            else
                z = z + zold3;
            end
        end
        if S2 == 1
            zold2 = z; 
        else
             z = z+ zold2;
        end
    end
    if S1 == 1
        zold1 = z; 
    else
        z = z+ zold1;
    end
end

But this way for performing the result is clearly verbose because there is a clear pattern from the third spin to N-1 spin, so what can I do ? I have been thinking different method for a long time like applying for to a function but couldn´t find any solution. Plz help. (P.D: The analitical solution is given by solving the 1D ising model for $\beta J=1$ with free boundary condition employing brackets method and is: $$Z = 2^N\cosh(1)^{N-1}$$ )