How can I draw a bifurcation diagram for the given code?

clc

close all

k1 = 3;

k_1 = 1;

k2 = 2.5;

k3 = 1;

k_3 = 1;

k4 = 2;

k5 = 1;

E1 = 1;

E2 = 2;

K1 = (k_1+k2)/k1;

K2 = (k_3+k4)/k3;

a = -(k5*E1)^2;

b = -((k3*E2+(k5*E1)^2)*(k1*(K1+E1)+k3*K2));

c= -(k3*E2*(k1*(K1+E1)+k4)+k1*(k5*E1)^2*(k2*E1+k3*K2*(K1+E1)));

d = -(k1*k3*(k4*K1*E2+E1*E2*(k2+k4)+k2*k5^2*K2*E1^3));

e=-(k5*E1);

f=-(k5*E1*(k1*K1+k3*(K2+E2)));

g = -(k1*k3*k5*E1*K1*(K2+E2));

m1 = -5*(a*e)^2;

m2 = (8*a^2*e*g-4*a^2*f^2-4*b^2*e^2+8*a*e^2*c);

m3 = (2*e*g-f^2)*(3*b^2-6*a*c)+3*e^2*(2*b*d-c^2)-3*a^2*g^2;

m4 = 2*(2*e*g-f^2)*(c^2-2*b*d)+2*g^2*(2*a*c-b^2)-2*e^2*d^2;

m5 = d^2*(2*e*g-f^2)+g^2*(2*b*d-c^2);

w = 0.001:0.01:0.9;

T = (m1*w^8+m2*w^6+m3*w^4+m4*w^2+m5)*((g-e*w^2)^2+(f*w)^2)*((d-b*w)^2+(c*w-a*w^3)^2);

Answer 1

[Chunru]

When do array operation, you should use .^ or .*:

k1 = 3;
k_1 = 1;
k2 = 2.5;
k3 = 1;
k_3 = 1;
k4 = 2;
k5 = 1;
E1 = 1;
E2 = 2;
K1 = (k_1+k2)/k1;
K2 = (k_3+k4)/k3;
a = -(k5*E1)^2;
b = -((k3*E2+(k5*E1)^2)*(k1*(K1+E1)+k3*K2));
c= -(k3*E2*(k1*(K1+E1)+k4)+k1*(k5*E1)^2*(k2*E1+k3*K2*(K1+E1)));
d = -(k1*k3*(k4*K1*E2+E1*E2*(k2+k4)+k2*k5^2*K2*E1^3));
e=-(k5*E1);
f=-(k5*E1*(k1*K1+k3*(K2+E2)));
g = -(k1*k3*k5*E1*K1*(K2+E2));
m1 = -5*(a*e)^2;
m2 = (8*a^2*e*g-4*a^2*f^2-4*b^2*e^2+8*a*e^2*c);
m3 = (2*e*g-f^2)*(3*b^2-6*a*c)+3*e^2*(2*b*d-c^2)-3*a^2*g^2;
m4 = 2*(2*e*g-f^2)*(c^2-2*b*d)+2*g^2*(2*a*c-b^2)-2*e^2*d^2;
m5 =  d^2*(2*e*g-f^2)+g^2*(2*b*d-c^2);
w = 0.001:0.01:0.9;
T = (m1*w.^8+m2*w.^6+m3*w.^4+m4*w.^2+m5).*((g-e*w.^2).^2+(f*w).^2).*((d-b*w).^2+(c*w-a*w.^3).^2);
plot(w, T)

[Akhtar Jan]

Thank you