Let $A = (Q = \{q_1, \dots, q_n\}, \Sigma, \delta, Q_F)$ be a (nondeterministic) finite automation with starting state $q_1$, $Q_F \subseteq Q$ and $\delta \subseteq Q\times\Sigma\times Q$.
Let $Q_i(z)$ the generating function for all the words that can be accepted starting in $q_i$, that is the $n$th coefficient of its series expansion $[z^n]Q_i = |\{w \mid |w| = n \wedge w \text{ accepted from } q_i\}|$.
Clearly:
$Q_i(z) = \left[ q_i \in Q_F \right] + \sum\limits_{(q_i, a, q_j) \in \delta} x \cdot Q_j(z)$
Solve the resulting (linear) equation system for $Q_1$ (using Mathematica or a similar tool). Then, $[z^n]Q_1$ is the desired quantity.
This goes back to a technique introduced for grammars by Chomsky and Schützenberger (1963); it easily transfers to finite automata.
Edit: If you want to account for $\varepsilon$-transitions, just leave out factor $x$ in the sum for the corresponding transition. Similiarly, if you have "compressed" edges, i.e. instead of symbol $a \in \Sigma$ a word $w \in \Sigma^k$ on a transition, replace $x$ with $x^k$.