function y=eigpowers(P,n) 
%EIGPOWERS express the n-th power of a 
% diagonalizable matrix as a sum of its eigenvalues.
% For example, given P with eignevales lambda_0, ...
%P^n = lambda0^n*A0 + lambda1^n*A1+ ... 
m = size(P,1); %determine the number of rows, or columns. 
[X,D] = eig(P); % determine eigenvectors and eigenvalues 
X_inv = inv(X); 
y(:,:,1) = D^n; % first matrix is D^n
for i=1:m, 
   sel = zeros(m); 
   sel(i,i) = 1; % sel is a matrix to deterine A_i 
   y(:,:,i+1) = X*sel*X_inv; 
end 

