Let $f(x) = x(x-2)$. When $f$ is applied to itself, all points that have a y of 3 stay at 3 since $3(3-2) = 3$. The function is centred at 1 so adding 1 to the x centers it at 0. This lead to the function, $g(x) = f^{n}(x+1)$, where the superscript is the number of iterations of the function onto itself. When the graph is stretched along the x-axis by a factor of $2^{n-1}$, lowered by 1, and then divided by 2, the function seems to approximate cosine. Why does this approximate cosine?
$$g(x) = \frac{(f^n(\frac{x}{2^{n-1}} + 1) - 1)}{2} \approx \cos(x)$$
Example, n = 10, the function on top is $g(x)$, on bottom is $cos(x)$:
The actual solution that I found to this problem actually didn't use the fact that $\cos(2x) = 2\cos(x)^2 - 1$, which I thought it would, but the solution is still quite fun nonetheless.
I started off by noticing that $f(x+1)-1$ is the same as $f(x)$ with respect to how it changes, but the computations are much easier. This should make sense since translations don't actually qualitatively affect a graph. I will now redefine $f(x)$ as $f(x+1)-1$ since I won't use the original throughout the rest of this proof. This means the new function for $f$ is $f(x) = (x+1)(x-1) - 1 = x^2-1-1 = x^2-2$.
Now, I will change the definition of $f$ even further by removing the need to divide by $2^{n-1}$ at the end. I will do this by saying that every time $f$ is applied, all the values of $x$ are mapped like so, $x \rightarrow \frac {x}{2}$.
The following is an example of how the function now works, $$f^1(x) = x^2 - 2$$
$$f^2(x) = \left( \left( \frac {x}{2} \right)^2 - 2 \right)^2 - 2 = \left( \frac {x^4} {16} - 4\left(\frac{x^2}{4}\right)+4\right)-2 = \frac{x^4}{16}-x^2+2$$
$$f^3(x)=\left( \frac{(\frac{x}{2})^4}{16}-\left(\frac{x}{2}\right)^2+2 \right)^2-2=\frac{x^8}{65536}-\frac{x^6}{512}+\frac{5x^4}{64}-x^2+2$$ $$\vdots$$
As you can see, for each iteration, the $x$ values of the function are all halved and the function is then applied. You may also have noticed that the coefficient of $x^2$ and the number $2$ stayed constant from the second iteration onwards. This is the main point of the proof, finding which number the coefficients converge to.
It is simple to prove that the constant always stays as $2$ after the first iteration. Since the constant has a power of $0$, the only thing that will multiply to give the constant is the constant itself. That means that the constant will always be $2(2)-2$, the minus two comes from the function, $f(x) = x^2 - 2$.
It is similarly simple to prove that the coefficient of $x^2$ always stays as $-1$, this is done in much the same way as for the constant. The things that multiply to give an $x^2$ term are the constant and $x^2$, that means that the $x^2$ term will be $2(a)\left(\frac{x}{2}\right)^2 + 2(a)\left(\frac{x}{2}\right)^2 = 4a\left(\frac{x}{2}\right)^2 = ax^2$, where $a$ is the previous coefficient. This means that when the constant becomes $2$, the coefficient of $x^2$ won't change, and when this happens, the coefficient of $x^2$ is $-1$.
What is more interesting is predicting what the coefficients of higher power terms will be. As an example, I will pick the next term, $x^4$. This time, the things that multiply to give an $x^4$ term are the $x^4$ term and the constant, and the $x^2$ term with itself. Again using $a$ as the previous coefficient of the $x^4$ term, the term has the following equation, $4a\left( \frac {x}{2} \right)^4 + (-1)(-1)\left( \frac {x}{2} \right)^4$. It will become obvious that we can continuously iterate this, inputting $a$ into itself to show that this does not depend on the original value of $a$ and can be iterated until all the coefficients before it have settled, meaning that their variability does not need to be taken into account.
The value for $a$ is represented by this expression, $\frac{a}{4} + \frac{1}{16} = \frac {4a + 1}{16}$. If this is iterated infinitely many times, inputting $a$ into itself, the following infinite fraction is produced,
$$\frac {4\left(\frac {4\left(\frac {4\left(\frac {4 \cdots + 1}{16}\right) + 1}{16}\right) + 1}{16}\right) + 1}{16}$$
A nicer way to represent this infinite fraction is the following, $a = \frac {4a + 1}{16}$ which encapsulates the idea of infinitely applying this to itself. This can also be solved quite simply, $a = \frac{4a + 1}{16}$ $\Rightarrow$ $16a = 4a + 1$ $\Rightarrow$ $12a = 1$ $\Rightarrow$ $a = \frac{1}{12}$. This means that the value of the coefficient for $x^4$ converges to $\frac{1}{12}$.
I will do one more example. For the $x^6$ term, the coefficient has the following equation which can be derived from its infinite fraction, $a = \frac{a - \frac{1}{6}}{64}$ which gives a value of $\frac {-1}{360}$ for the coefficient of the $x^6$ term.
We now need a way to calculate the value of all the terms instead of going through them one by one. The coefficient will always be the sum of all the products of terms where the sum of their two degrees adds to the desired degree. For any degree $2n$, these terms will be $\frac{1}{2^{2n}}\sum_{k=1}^{n-1} \frac{2}{(2k)!} \frac{2}{(2n-2k)!} = \frac{1}{2^{2n}}\sum_{k=1}^{n-1} \frac{4}{(2k)!(2n-2k)!}$. The $\frac{1}{2^{2n}}$ is there because the $x$ values are getting halved and the degree, when multiplied together, is $2n$. The first and last terms aren't included since they contain $a$, this amounts to $\frac{4a}{2^{2n}}$. Finally, we have to take into account that the sign of the coefficients flips every term. This gives the final equation,
$$a = \frac{4a}{2^{2n}} - \frac{(-1)^n}{2^{2n}}\left( \sum_{k=1}^{n-1} \frac{4}{(2k)!(2n-2k)!} \right)$$
$$(2^{2n}-4)a = - (-1)^n\left( \sum_{k=1}^{n-1} \frac{4}{(2k)!(2n-2k)!} \right)$$
$$(2^{2n}-4)a = - (-1)^n\left(\frac{-8}{(2n)!} + \sum_{k=0}^{n} \frac{4}{(2k)!(2n-2k)!} \right)$$
Note: I did not assume the answer in the question, the factorial representation builds upon itself since it holds for the first and second and this shows that it forces it to hold for all the others.
Solving the sum in the equation is pretty simple, since, $$(x + 1)^n = \sum_{k=0}^{n} \frac{n!x^n}{(k)!(n-k)!}$$ $$(1 + 1)^{2n} + (-1 + 1)^{2n} = 2^{2n} + 0 = 2^{2n}$$ $$=\sum_{k=0}^{2n} \frac{(2n)!1^k}{(k)!(2n-k)!} + \sum_{k=0}^{2n} \frac{(2n)!(-1)^k}{(k)!(2n-k)!} = \sum_{k=0}^{n} \frac{2(2n)!}{(2k)!(2n-2k)!}$$ $$2^{2n+1} = \sum_{k=0}^{n} \frac{4(2n)!}{(2k)!(2n-2k)!}$$ $$\frac{2^{2n+1}}{(2n)!} = \sum_{k=0}^{n} \frac{4}{(2k)!(2n-2k)!}$$
This can now be substituted back into the previous equation,
$$(2^{2n}-4)a = - (-1)^n\left(\frac{-8}{(2n)!} + \frac{2^{2n+1}}{(2n)!} \right)$$ $$a = \frac{-(-1)^n\left(-8 + 2^{2n+1} \right)}{(2n)!(2^{2n}-4)}$$ $$a = \frac{2(-1)^n\left(2^{2n} - 4\right)}{(2n)!(2^{2n}-4)}$$ $$a = \frac{2(-1)^n}{(2n)!}$$
Finally, by dividing by two, which is the only step that was not included in the new function, the cosine function is produced.