One way to do this is using Jordan normal form. It works well if you can write down the eigenvectors of $A$ explicitly. But there is another way. Consider the matrix generating function
$(I - At)^{-1} = \sum_{n \ge 0} A^n t^n.$
The coefficient of $t^n$ is precisely a matrix whose values describe the number of paths of length $n$ between the vertices of $G$. On the other hand, this is just the inverse of a matrix $I - At$, so its entries can be computed using the adjugate, which gives
$(I - At)^{-1} = \frac{\text{adj}(I - At)}{\det(I - At)}$
hence
$\sum_{n \ge 0} (A^n)_{ij} t^n = \frac{ \text{adj}(I - At)_{ij} }{\det(I - At)}.$
But $\text{adj}(I - At)_{ij}$ is easy to calculate. Thus the generating function for the sequence you want can be explicitly written as the ratio of two polynomials which are easy to compute given $A$, and then you can use partial fraction decomposition to turn this into an explicit formula.
If you are less interested in an explicit formula than just efficiently computing $(A^n)_{ij}$, then it may be more productive to just compute $A^n$ using binary exponentiation.