%% This code produces one-parameter bifurcation diagram for the 3D Lactose Operon Model
function Code_for_Ex_2_4
clc
ssA       = [10^(-10):10^(-3):1e0];
init      = [1 1 1]*0.006;
for i=1:length(ssA)
    As              = ssA(i);
    SSValues(i,:)   = fsolve(@bifur,init,[],As);    
    ssM             = SSValues(:,1);
    ssB             = SSValues(:,2);
    ssL             = SSValues(:,3);
    init            = SSValues(i,:) ;
end
%Below Plot ssL versus ssA, choose y-axis in log scale
% .........................
axis([0.001 0.7 0.001 1])


%% Steady State Equations
function dy=bifur(y,As)


%% Parameters
% Enter parameter values below from the 3D lac model paper published in Chaos
% ...............

dA = alphaA*B*L/(KL+L)-betaA*B*A/(KA+A) - (gammaA+mu)*A;
dM = alphaM*((1+K1*exp(-n*tauM*mu)*A^n)/(K+K1*exp(-n*tauM*mu)*A^n))  - (gammaM+mu)*M;
dB = alphaB*exp(-tauB*mu)*M -(gammaB+mu)*B;
dy = [dA dM dB]';

