setdefault

% set k here
k = 1980;

llim = 1.5;
lambda = linspace( -llim, llim, 300 );
del = zeros( 1, length(lambda) );

for j = 1:length(lambda)
    % First find the steady state
    M = [   r/k    r/k       b;
            r/k*a  r/k       b-eps;
            c*b    c*(b-eps) 0       ];
    ss = inv(M) * [r;r;d];
    
    Ja = [   -r/k*ss(1)-lambda(j)   -r/k*ss(1)                -b*ss(1);
             -r/k*a*ss(2)           -r/k*ss(2)-lambda(j)      -(b-eps)*ss(2);
             c*b*ss(3)               c*(b-eps)*ss(3)          -lambda(j)           ];
    del(j) = det(Ja);
end

close all;
plot(lambda,del);
xlabel('lambda');ylabel('det(J-lambda*I)');
title('Characteristic equation of the linearization');
grid on;

% Finding the root numerically
% We assume that the characteristic function is initially positive
for j = 2:length(lambda)
    if(del(j)<0)
        j_crit = j;
        break;
    end
end

% Report the results of search
crit_m1 = [lambda(j_crit-1) del(j_crit-1)]
crit_m0 = [lambda(j_crit)   del(j_crit)  ]
avg = 1/2*(lambda(j_crit-1)+lambda(j_crit))

% Linearized interpolation 
avg = lambda(j_crit) + -del(j_crit)*(lambda(j_crit)-lambda(j_crit-1))/(del(j_crit)-del(j_crit-1))