The orientation-preserving symmetry group $G_n$ of the $n$-dimensional cube is an index two subgroup of the full symmetry group, which is $S_n\times\{\pm 1\}^n$. By the Polya-Burnside counting theorem, the number of ways to color the cube with two colors is just the average size of $2^{C}$, where $C$ is the number of cycles in the action of a random element $g$ of $G_n$ on the faces of the cube. Write $g = (\sigma, (s_i)_{i=1,...,n})$, where $\sigma$ is a permutation of the basis vectors and $(s_1,...,s_n)$ is a collection of signs such that $\prod_{i\le n} s_i$ is the sign of $\sigma$.
The strategy is to now fix the permutation $\sigma$ and average over collections of signs $(s_1, ..., s_n)$. Let $c$ be the number of cycles of $\sigma$, so that the sign of $\sigma$ is $(-1)^{n-c}$. A cycle $(i_1, ..., i_k)$ of $\sigma$ either stays a cycle of $g$ or splits into two cycles of $g$ depending on whether an odd or even number of the $s_{i_k}$'s are equal to $-1$. If we set $p_m$ to be the product of the $s_{i_k}$'s in the $m$th cycle, we see that $2^C = \prod_{m\le c} (3+p_m)$, and $\prod_{m\le c}p_m = (-1)^{n-c}$. The average of $\prod_{m\le c} (3+p_m)$ over all choices of $p_1, ..., p_c$ satisfying $\prod_{m\le c}p_m = (-1)^{n-c}$ comes out to $3^c + (-1)^{n-c}$ (easy exercise).
Since the average of $(-1)^{n-c}$ is $0$ for $n>1$ (i.e. exactly half of all permutations are odd), we just need to find the average value of $3^c$ where $c$ is the number of cycles of a random permutation $\sigma$ of $S_n$. By Polya-Burnside again, the expected size of $3^c$ is exactly the number of ways to color a set of size $n$ with $3$ colors, up to permutation of the $n$ elements, and this is just $\binom{n+3-1}{3-1} = \frac{(n+1)(n+2)}{2}$.
In general, the number of ways to paint an $n$-dimensional cube with $k$ colors comes out to
$\binom{n+\frac{k^2+k}{2}-1}{\frac{k^2+k}{2}-1}+(-1)^n\binom{n-\frac{k^2-k}{2}-1}{-\frac{k^2-k}{2}-1}.$
The second summand can be omitted if you allow orientation-reversing symmetries of the cube.
