This is matlab function to calculate sum of A to N power using matrix diagonalization, it assumes that matrix is a square matrix.
function [x] = powersum(A, m)
% Greg Bugaj
% Y = powersum(A, n) gives an sum of matrices to N power, if matrix size
% is less than 2 then we simply return the input matrix
% Input : A - an nxn matrix
% n - How many matrices to sum
%Output x - summed matrix to N power.
% P is our eigenvector, d is the diagonal to verify( inv(p)*A*p )
x = A;
[p d] = eigs(A);
for i = 2 : m
% x = x+(A^i); -- Just raise to power and sum
% note b * inv(A) is same as b/p
x = x + ((p * d^i) / p);
end