The general expression of the Taylor expansion for functions with two variables is
f(x,h)=sum_(n=0)^oo(1/(n!) sum_(k=0)^n((n),(k))((partial^nf)/(partialx^(n-k)partialh^k))_({x_0,h_0})(x-x_0)^(n-k)(h-h_0)^k)
In this case, considering x_0=0,h_0=0
((partial^nf)/(partialx^(n-k)partialh^k))_(x=0,h=0) = cos((n+k)pi/2)
so
cos(x+h) = sum_(n=0)^oo(1/(n!) sum_(k=0)^n((n),(k))cos((n+k)pi/2)x^(n-k)h^k)
but
cos((n+k)pi/2) = i^(n+k)((1+(-1)^(n+k))/2) and
((n),(k)) = (n!)/((n-k)!k!)
so
cos(x+h) = sum_(n=0)^oosum_(k=0)^n ( i^(n+k)((1+(-1)^(n+k))/2))/((n-k)!k!)x^(n-k)h^k
with i = sqrt(-1)
Another way to do that is knowing that from
cosx = sum_(k=0)^oo (-1)^k(x^(2k))/(2k!) follows
cos(x+h) = sum_(k=0)^oo(-1)^k((x+h)^(2k))/(2k!) but here the variables x+h appear added.
Multivariate series handling is very cumbersome because the required notation effort needed.