%
%  math495/stat490 -- 03 nov 03 -- djm
%
%  w10ehren.m:  

PT = [0 1 0 0 0 ;
      4 0 2 0 0 ; 
      0 3 0 3 0 ;
      0 0 2 0 4 ;
      0 0 0 1 0 ] /4;
 
P = transpose(PT);

%  eigensolution
[evec,eval] = eig(PT);

%  basis matrix
v1 = [1/16 1/4 3/8 1/4 1/16]';
v2 = [1 -4 6 -4 1]';
v3 = [1 0 -2 0 1]';
v4 = [1 2 0 -2 -1]';
v5 = [-1 2 0 -2 1]';

V = [v1 v2 v3 v4 v5];

%  example
y0 = [1 0 0 0 0]';

a = V\y0;