% ex4

op=optimset('TolX',1e-10,'TolFun',1e-8,'LargeScale','on');
p=0.1309;
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);

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

dx1 = [0.8 1];
dx2 = [3 6];
figure(1);
axis([dx1 dx2]);
hold on
c='gr';

op=odeset('AbsTol',1e-10,'RelTol',1e-8);

for i=1:50
  x0 = rand(2,1).*([dx1(2);dx2(2)]-[dx1(1);dx2(1)])+[dx1(1);dx2(1)];
  [t,x]=ode15s('ex4', [0, 20], x0, op, p);
  plot(x(:,1),x(:,2));
end

if isreal(L1)
  [z,i]=sort(abs(diag(L1)));
  plot([-X1(1,1) X1(1,1)]+xs1(1),[-X1(2,1) X1(2,1)]+xs1(2),c(i(1)),[-X1(1,2) X1(1,2)]+xs1(1),[-X1(2,2) X1(2,2)]+xs1(2),c(i(2)));
end

if ns > 1
  if isreal(L2)
    [z,i]=sort(abs(diag(L2)));
    plot([-X2(1,1) X2(1,1)]+xs2(1),[-X2(2,1) X2(2,1)]+xs2(2),c(i(1)),[-X2(1,2) X2(1,2)]+xs2(1),[-X2(2,2) X2(2,2)]+xs2(2),c(i(2)));
  end
  plot(xs2(1),xs2(2),'ro');
end

if ns > 2
  if isreal(L3)
    [z,i]=sort(abs(diag(L3)));
    plot([-X3(1,1) X3(1,1)]+xs3(1),[-X3(2,1) X3(2,1)]+xs3(2),c(i(1)),[-X3(1,2) X3(1,2)]+xs3(1),[-X3(2,2) X3(2,2)]+xs3(2),c(i(2)));
  end
  plot(xs3(1),xs3(2),'ro'); 
end

plot(xs1(1),xs1(2),'ro');
plot(dx1,[0 0],'k',[0 0],dx2,'k');
xlabel('x_1'); ylabel('x_2');
