From Sage v.8.0.beta3,
sage: A = Matrix([[1,2],[3,4]])
sage: n = SR.var('n')
sage: A^n
[ 1/22*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 11) - 1/22*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) - 11) 2/33*sqrt(33)*(1/2*sqrt(33) + 5/2)^n - 2/33*sqrt(33)*(-1/2*sqrt(33) + 5/2)^n]
[-1/88*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 11)*(sqrt(33) - 3) - 1/88*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 3)*(sqrt(33) - 11) 1/66*sqrt(33)*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 3) + 1/66*sqrt(33)*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) - 3)]
which is
(122(−12√33+52)n(√33+11)−122(12√33+52)n(√33−11)233√33(12√33+52)n−233√33(−12√33+52)n−188(−12√33+52)n(√33+11)(√33−3)−188(12√33+52)n(√33+3)(√33−11)166√33(12√33+52)n(√33+3)+166√33(−12√33+52)n(√33−3))
Also in the case of 2x2 matrices there is a derivation of the closed formula for An in this note by K. S. Williams.