Self-defined zeta function never terminates

Question

I am trying to build a program to compare the partial sum of the Riemann Zeta function to the built-in Matlab function zeta(s). I want the function to output the minimum number of terms to achieve an accuracy of 0.1 percent. I thought a while loop would be the best approach, but my program is running so slow; I have yet to get a result from it.

function[n] = riemannzeta(s)
error = 1; n = 1; an = 1; S = an;
while error >= 0.1
    an = 1/n^s;
    S = S + an;
    n = n + 1;
    z = zeta(s); 
    error = ((S - z)/z)*100;
end
end

I call it with:

riemannzeta(3)

Answer 1

The main problem is that your definition of the zeta function is wrong because you initialize the value of the sum to 1 and then add 1 in the first step. You either need to initialize at 0 or start the loop at 1/2^s. You also need to take the absolute value of the error.

Here is the start-at-two version:

function n = riemannzeta(s)
error = 1; n = 1; an = 1; S = 1;
z = zeta(s); 
while error >= 0.001
    n = n + 1;
    an = 1/n^s;
    S = S + an;
    error = abs(S - z)/z;
end
end

If I run riemannzeta(3) I get a value of 20.