% ex4
function L=rloc_ex4(p)

op=optimset('TolX',1e-10,'TolFun',1e-8,'LargeScale','on');
xs=fsolve('ex4_ss',0+(p>0.106),op,p);

if p == 0.0889 | p == 0.10574
   ns = 2;
else if p > 0.0889 & p < 0.10574
      ns = 3;
   else
      ns = 1;
   end
end
      
xs1 = [xs; 14/3*xs];
A = [-1-p*exp(xs1(2)) p*(1-xs1(1))*exp(xs1(2));
     -14*p*exp(xs1(2)) -3+14*p*(1-xs1(1))*exp(xs1(2))];
[X1,L1] = eig(A);

L2=diag([0,0]);
L3=L2;

if ns > 1,
  xs=fsolve('ex4_ss',0.5+5*(ns<3),op,p);  
  xs2 = [xs; 14/3*xs];
  A = [-1-p*exp(xs2(2)) p*(1-xs2(1))*exp(xs2(2));
     -14*p*exp(xs2(2)) -3+14*p*(1-xs2(1))*exp(xs2(2))];
  [X2,L2] = eig(A);
end

if ns > 2,
  xs=fsolve('ex4_ss',5,op,p);  
  xs3 = [xs; 14/3*xs];
  A = [-1-p*exp(xs3(2)) p*(1-xs3(1))*exp(xs3(2));
     -14*p*exp(xs3(2)) -3+14*p*(1-xs3(1))*exp(xs3(2))];
  [X3,L3] = eig(A);
end

L=[diag(L1), diag(L2), diag(L3)];
