How to implement Steffenson method using Goertzel method for evaluating function values? The task is to implement function in Matlab that will solve numerically following function: 
$\sum _{k=0}^na_k\cos \left(kx\right)$
 using Steffenson method and Goertzel for evaluating the series. 
I have working Steffenson method part, but no idea how to impelement Goertzel and how to combine it into.
I'll be glad for any hint. 
 A: The Goertzel method is a generalization of the Horner scheme to evaluate trigonometric polynomials, originally for the extraction of single or few amplitudes from a Fourier transform. It essentially uses the trigonometric identity $$ \cos((k+1)x)+\cos((k-1)x)=2\cos x\cos(kx)$$
Inserting this identity into the expression gives
\begin{align}
\sum_{k=m}^{n}a_{k}\cos(kx)
&=a_0+a_1\cos(x)+\sum_{k=2}^na_k[2\cos x\cos((k-1)x)-\cos((k-2)x)]
\\
&=a_0-a_1\cos(x)+2\cos(x)\sum_{k=0}^{n-1}a_{k+1}\cos(kx)-\sum_{k=0}^{n-2}a_{k+2}\cos(kx).
\end{align}
Thus setting $$s_m(x)=\sum_{k=0}^{n-m}a_{k+m}\cos(kx)$$
we get $$s_0(x)=a_0-a_1\cos(x)+2\cos(x)s_1(x)-s_2(x)$$ for the above equation and in general
$$
s_{m-1}(x)=a_{m-1}+\cos(x)(2s_m(x)-a_m)-s_{m+1}(x)
$$
which gives the implementation
def Goertzel(a,x):
    n = len(a)-1
    cx = cos(x)
    s1, s0 = a[n], a[n-1]+a[n]*cx
    for k in range(n-2,-1,-1): # n-2 down to 0
        s1, s0 = s0, a[k] + (2*s0-a[k+1])*cx - s1
    return s0

If the compiler recognizes the multiplication with 2 and implements it as a simple increment in the binary exponent, the cost of this algorithm is the evaluation of the cosine, once, and $n$ multiplication, in total not much more than a polynomial evaluation of the same degree. 
