You can expand() to get a 'longer' form
syms a b;
f = cos(a) + cos(b);
expand(f^2)
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.
. 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]))
test2 = eval(subs(g, [a b],[1 2]))
