clear;
beta=0.9;
alpha=0.75;
sigma=1;
delta=0.3;
rho=0.95;

kbar=((1-(1-delta)*beta)/(alpha*beta))^(1/(alpha-1));
cbar=kbar^alpha-delta*kbar;
ybar=kbar^alpha;
A(1,1)=1-(1-delta)*beta;
A(1,2)=(1-(1-delta)*beta)*(alpha-1);
A(1,3)=-sigma;
A(2,2)=kbar;
A(3,1)=1;
B(1,3)=-sigma;
B(2,1)=ybar;
B(2,2)=alpha*ybar+(1-delta)*kbar;
B(2,3)=-cbar;
B(3,1)=rho;

C=inv(B)*A;

[ve,MU]=eig(C);                %ve is a matrix with eigen vectors as columns and MU is a matrix eigen values on the diagonal and zeros elsewhere

if (MU(1,1)<1 & MU(2,2)<1) | (MU(1,1)<1 & MU(3,3)<1) | (MU(2,2)<1 & MU(3,3)<1) | (MU(1,1)>1 & MU(2,2)>1 & MU(3,3)>1)
    display('No saddle path stability.');
else
    
        val= diag(MU);
        [val ve'];
        t=flipud(sortrows([val ve']));  
        MU=diag(t(:,1));
        ve=t(:,2:4);
        P=inv(ve')      
        
        MU1=MU(1:2,1:2);
        MU2=MU(3,3);
        P11=P(1:2,1:2);
        P12=P(1:2,3);
        P21=P(3,1:2);
        P22=P(3,3);
        
        -inv(P22)*P21
        inv(P11-P12*inv(P22)*P21)*inv(MU1)*(P11-P12*inv(P22)*P21)
    end