Hi ELISHA ANEBI,
I noticed that in this equation 'H1=C2*c5*c6*c8*c9;' the C2 sholud be changed to c2. Additionally, there is a syntax error in the line for t=1:1:len;. The semicolon at the end should be removed.
To create the bifurcation diagram, you can modify the code as follows:
Rev_value=0.018:0.01:4;
Root_array=zeros (length (Rev_value), 2);
qI=0.001923; qA=0.00000004013; etaA=0.1213; etaQ=0.003808; w=0.5925;
Lambda=0.1598643e-7; theta=0.022;
delta=0.125; mu=0.01119; pi=0.464360344; deltaQ=6.847e-4; beta=1.086e-1;
qE=1.8113e-4; rhoQ=0.0815; a=0.16255; k=0.15; v1=0.71; v2=0.29;alpha=0.57e-1; deltaI=0.00000000223; rhoA=0.1; rhoI=0.0666666;
c1=mu+v1; c2=1-w; c3=mu+alpha+v2;c4=1-k; c5=qE+delta+mu; c6=rhoA+mu; c7=delta*(1-a); c8=rhoI+qI+deltaI+mu; c9=rhoA+deltaQ+mu;
b1=1-theta;
B=delta*a*c8*c9+c7*k*qI*rhoQ+c8*k*qE*rhoQ;
G=qI*c7*c6+qE*c8*c6;
H1=c2*c5*c6*c8*c9;
H2=c5*c6*c8*c9*(c3*+c1*c2)-(b1*c2+c2*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
H3=c5*c6*c8*c9*(c3*c1-v1*v2)*(1-b1*Rev_value)-(c3*theta+c2*v1*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
hold on
for i = 1:length(Rev_value)
Rev = Rev_value(i);
% Calculate G
G = qI * c7 * c6 + qE * c8 * c6;
% Coefficients of quadratic equation H1, H2, H3
H1 = c2 * c5 * c6 * c8 * c9;
H2 = c5 * c6 * c8 * c9 * (c3 + c1 * c2) - (b1 * c2 + c2 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
H3 = c5 * c6 * c8 * c9 * (c3 * c1 - v1 * v2) * (1 - b1 * Rev) - (c3 * theta + c2 * v1 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
% Coefficients of quadratic equation p = [H1, H2, H3]
p = [H1, H2, H3];
r = roots(p);
len = length(r);
for t = 1:len
% Store the real part of the roots in the Root_array
if isreal(r(t))
Root_array(i, t) = real(r(t));
end
end
end
% Plot the bifurcation diagram
plot(Rev_value, Root_array, '.')
xlabel('Rev')
ylabel('Roots')
title('Bifurcation Diagram')
The code will generate a bifurcation diagram by plotting the real part of the roots against the parameter Rev.
