how to expand/rewrite/simplify (cos(a)+cos(b))^2?

I have the code below. How to expand/rewrite/simplify the expression f^2 into cos(2a), cos(2b), cos(a+b), cos(2a), cos(2b), cos(a-b), cos(b-a)?

==========================================

syms a b;
f=cos(a)+cos(b);
combine(f^2)

ans = 

shows (cos(a) + cos(b))^2

rewrite(f^2, 'sincos') 

ans = 

shows (cos(a) + cos(b))^2

simplify(f^2)

ans = 

shows (cos(a) + cos(b))^2

==========================================

[Nathan Hardenberg]

You can expand() to get a 'longer' form

syms a b;
f = cos(a) + cos(b);
expand(f^2)

ans = 

I don't know exactly how you want to rewite it. From your question I understand that you want to get

But simply testing numerically you can see that . So can you please elaborate on what you want your result to be, since I most likely misunderstood it.

g = cos(2*a) * cos(2*b) * cos(a+b) * cos(2*a) * cos(2*b) * cos(a-b) * cos(b-a);
test1 = eval(subs(f^2, [a b],[1 2]))

test1 = 0.0154

test2 = eval(subs(g, [a b],[1 2]))

test2 = -0.0214

[cfy30]

The result I am looking for is without the square, the last expression below. Any way to do it?

(cos(a) + cos(b))^2

=cos(a)^2 + 2cos(a)cos(b) + cos(b)^2

=1/2*[1+cos(2a)] + cos(a+b)+cos(a-b) + 1/2*[1+cos(2b)]

[Nathan Hardenberg]

I could be wrong, but I don't think this is possible. In the rewrite() documentation you can find a list of functions that can be replaced, but squares/powers are not (yet?) one of the options.

Still I want to mention that you can find multiple simplifications with the following Code:

syms a b
expr = cos(a) + cos(b);
simplify(expr,'Steps',15,'All',true)

ans =