

close all;



%parameters

a = 1.5;

b = 0.01;

c = 0.5;

d = 1.0;

r = 1.0;

eps = 0.009;



%fix k here

k = 400;



M = [r/k r/k b; r/k*a r/k (b-eps); c*b c*(b-eps) 0];

% fp: fixed point solution

fp = inv(M)*[r;r;d]



lim = 2.5;

lambda = linspace(-lim,lim,600);

de = zeros(1,length(lambda));



for j=1:length(lambda)

   M = [-r/k*fp(1)-lambda(j) -r/k*fp(1) -b*fp(1);

        -r/k*a*fp(2) -r/k*fp(2)-lambda(j) -(b-eps)*fp(2);

        c*b*fp(3) c*(b-eps)*fp(3) -lambda(j)];

   de(j) = det(M);

   plot(lambda(j),de(j),'.k');

   hold on;

end



grid on;

xlabel('lambda');

ylabel('det(J-lambda*I)');

axis([-lim lim -1 1]);



hold on;

plot(lambda,zeros(1,length(lambda)));



de