Monodromy matrix coding not working

% Initialization

clc,clear,cnt=1;

syms d a phia;

syms tau x0_1 x0_2 x0_3 x01 x0 x0_4;

x0=[x0_1;x0_2;x0_3;x0_4];

Vin=100;L1=200e-6;L2=200e-6;C=20e-6;R=57.6/2;T=10e-6;Ki=500;Kp=5;Kil=1/4;Kvc=5/240;mc=-0.3;

iL10=0;iL20=0;Vc0=0;Vp0=0;Vref1=5;a1=0;

A_on_on=[-1/R/C 0 0 0;0 0 0 0;0 0 0 0;Ki*Kvc 0 0 0];

A_on_off=[-1/R/C 0 1/C 0;0 0 0 0;-1/L2 0 0 0;Ki*Kvc 0 0 0];

A_off_on=[-1/R/C 1/C 0 0;-1/L1 0 0 0;0 0 0 0;Ki*Kvc 0 0 0];

A_off_off=[-1/R/C 1/C 1/C 0;-1/L1 0 0 0;-1/L2 0 0 0;Ki*Kvc 0 0 0];

Ba=[0 0 0 0;0 0 1/L1 0;0 0 1/L2 0;0 0 0 -Ki];Bb=Ba;Bc=Bb;Bd=Bc;

Ua=[0;0;Vin;Vref1*(1+a1*sin(4*pi*tau/T-4*pi*d))];

Ub=[0;0;Vin;Vref1];

Vin1=80:1:125;

hold on;

for ii=38:42

%*****************************************

Vin=Vin1(ii);

sim('Interleaved_close_loop_supervision_control_slope_compensation');

iL100=iL1n(end-2,2);

iL200=iL2n(end-2,2);

Vc00=Vc1(end-2,2);

int00=Int(end-2,2);

x00=[Vc00;iL100;iL200;int00];

%duty cycle calculation

duty=Dutycycle1(end-1,1)/T;

B1=Ba;U=Ua;

if duty<0.5

%%%%%%%%%%%%%%When d<0.5;

A1=A_on_off;A2=A_off_off;A3=A_off_on;A4=A_off_off;

phi1=expm(A1*d*T);phi2=expm(A2*0.5*(1-2*d)*T);

phi3=expm(A3*d*T);phi4=expm(A4*0.5*(1-2*d)*T);

I1=int(expm(A1*(d*T-tau))*B1*U,tau,0,d*T);

I2=int(expm(A2*(0.5*T-tau))*B1*U,tau,0.5*T,(0.5+d)*T);

I3=int(expm(A3*((0.5+d)*T-tau))*B1*U,tau,0.5*T,(0.5+d)*T);

I4=int(expm(A4*(T-tau))*B1*U,tau,(0.5+d)*T,T);

else if duty>0.5

%%%%%%%%%%%%%%%%%%%%%%When d>0.5

A1=A_on_on;A2=A_on_off;A3=A_on_on;A4=A_off_on;

phi1=expm(A1*0.5*(2*d-1)*T);phi2=expm(A2*(1-d)*T);

phi3=expm(A3*0.5*(2*d-1)*T);phi4=expm(A4*(1-d)*T);

I1=int(expm(A1*(0.5*(2*d-1)*T-tau))*B1*U,tau,0,0.5*(2*d-1)*T);

I2=int(expm(A2*(0.5*T-tau))*B1*U,tau,0.5*(2*d-1)*T,0.5*T);

I3=int(expm(A3*(d*T-tau))*B1*U,tau,0.5*T,d*T);

I4=int(expm(A4*(T-tau))*B1*U,tau,d*T,T);

else if duty==0.5

%%%%%%%%%%%%%%%%%%%%%%When d=0.5

%A1=A_on_off;A2=A_off_on;A3=A_on_off;A4=A_off_on;

end

end

end

%--------------------------------

x1=phi1*x0+I1;

x2=phi2*x1+I2;

x3=phi3*x2+I3;

x4=phi4*x3+I4;

x1=subs(x1,[d,x0_1,x0_2,x0_3,x0_4],[duty,Vc00,iL100,iL200,int00]);

x2=subs(x2,[d,x0_1,x0_2,x0_3,x0_4],[duty,Vc00,iL100,iL200,int00]);

x3=subs(x3,[d,x0_1,x0_2,x0_3,x0_4],[duty,Vc00,iL100,iL200,int00]);

x4=subs(x4,[d,x0_1,x0_2,x0_3,x0_4],[duty,Vc00,iL100,iL200,int00]);

phi1=subs(phi1,d,duty);phi2=subs(phi2,d,duty);

phi3=subs(phi3,d,duty);phi4=subs(phi4,d,duty);

if duty<0.5

%%%%%%%%%%%%%%%%%%%%%%When d<0.5;

spa=-Kp*Kvc*(x1(3)*R-x1(1))/R/C-Kil*Vin/L1+Ki*(Vref1-Kvc*x1(1));

spb=-Kp*Kvc*(x3(2)*R-x3(1))/R/C-Kil*Vin/L2+Ki(Vref1-Kvc*x3(1));

sa=mc/T;

S12=[1-Kp*Kvc*x1(2)/C/(spa+sa) -Kil*x1(2)/C/(spa+sa) 0 x1(2)/C/(spa+sa);

Kp*Kvc*x1(1)/L1/(spa+sa) 1+Kil*x1(1)/L1/(spa+sa) 0 -x1(1)/L1/(spa+sa);

0 0 1 0;

0 0 0 1];

S34=[1-Kp*Kvc*x3(3)/C/(spb+sa) 0 -Kil*x3(3)/C/(spb+sa) x3(3)/C/(spb+sa);

0 1 0 0;

Kp*Kvc*x3(1)/L2(spb+sa) 0 1+Kil*x3(1)/L2/(spb+sa)-x3(1)/L2/(spb+sa);

0 0 0 1;

M=phi2*S12*phi1*phi4*S34*phi3

eig(M);

else if duty>0.5

%%%%%%%%%%%%%%%%%%%5When d>0.5

spa=Kp*Kvc*x1(1)/R/C-Kil*Vin/L1+Ki*(Vref1-Kvc*x1(1));

spb=Kp*Kvc*x3(1)/R/C-Kil*Vin/L2+Ki*(Vref1-Kvc*x3(1));

sa=mc/T;

S12=[1-Kp*Kvc*x1(3)/C/(spb+sa) 0 -Kil*x1(3)/C/(spb+sa) x1(3)/C/(spb+sa);

0 1 0 0;

Kp*Kvc*x1(1)/L2/(spb+sa) 0 1+Kil*x1(1)/L2/(spb+sa)-x1(1)/L2/(spb+sa);

0 0 0 1];

S34=[1-Kp*Kvc*x3(2)/C/(spa+sa) -Kil*x3(2)/C/(spa+sa) 0 x3(2)/C/(spa+sa);

Kp*Kvc*x3(1)/L1/(spa+sa) 1+Kil*x3(1)/L1/(spa+sa) 0 -x3(1)/L1/(spa+sa);

0 0 1 0;

0 0 0 1;

M =phi2*S12*phi1*phi4*S34*phi3

eig(M);

else

if duty==0.5

%%%%%%%%%%%%%%%%%%%%When d=0.5

A1=A_on_off;A2=A_off_on;A3=A_on_off;A4=A_off_on;

end

end

end

plot(real(eig(M)),imag(eig(M)),'o','MarkerSize',10,'color','g');

end