[SOLVED] What is the best software to simulate a boost converter with current control?

/%% the Data: value of the Boost converter´s element
%% christophe basso 225
Vin=20; % voltage input
Vout=40; % voltage output
r_cf=0.1 ;     %Output capacitor ESR
r_lf =0.1 ;    %inductr ESR
R= 200;       %Load resistor
C= 10.7e-6;         %Output Capacitor
L= 20e-6;          %inductor
M= Vout/Vin;          %conversion ration Vout/Vin
D=1-1/M;          %on duty cycle
D_=1-D;           %Off duty cycle, also denoted by 1 – D, depending on the mood
Fsw=2.2e6;        % switching  frequency
Tsw=1/Fsw;        %Switching period
Ri=1;         %sense resistor in current mode
S1=Vin/L;         %(S_1)inductor on slope, Vin/L for a boost converter
S2=(Vout-Vin)/L;         %inductor off slope, (Vout-Vin)/L
Sa=0.5*S2;         %external ramp compensation slope
mc=1+Sa/S1;         %compensation ramp according to ridley notation [2]:mc=1+Se/Sn;
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% it must under a condition to awitch between the DCM and CCM if Icrit>IL or k_crit
%%Current mode, CCM(reference equations 1):
k=(R/(2*Ri))*(Vin/Vout); %% DC gain
b=(((1+(Sa/S1))*D_)-0.5)/Fsw;
a=1/((pi*Fsw)^2);
w_z1=1/(r_cf*C);
w_z2=(R*(D_^2))/L; % Righ half plane zero
w_p1=2/(R*C);

%G(s)=Vout(s)/Verr(s)
%G1
G1=tf(k*[1/w_z1 1], [1/w_p1 1]);
%G2
G2=tf([-1/w_z2 1],[a b 1]);

G_ccm=G1*G2;
bode(G1,{10,1e7});
bode(G2,{10,1e7});
options = bodeoptions;
options.FreqUnits = 'Hz'; % or 'rad/second', 'rpm', etc.
bode(G_ccm,options,{10,1e7})
grid on
%%
%%Vout(s)/Vin(s)
gain=Vout/(2*Vin);
G_in_o=tf(gain*[1/w_z1 1],[1/w_p1 1]);
bode(G_in_o,{10,1e7})
grid on
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%current mode, DCM (reference equations 1)
tau=L/(R*Tsw);
M=1+sqrt(1+(2*(D^2))/tau)/2;
Fm=1/(S1*mc*Tsw);
Hd= ((2*Vout)/D)*(M-1)/(2*M-1);

w_z1d=1/(r_cf*C);
w_z2d= R/((M^2)*L);
w_p1d=(2*M-1)/((M-1)*(R*C));
w_p2d=2*Fsw*((1-(1/M))/D)^2;
G1d=tf(Fm*Hd*[1/w_z1d 1],[1/w_p1d 1]);
G2d=tf([-1/w_z2d 1],[1/w_p2d 1]);
G_dcm=G1d*G2d;

options = bodeoptions;
options.FreqUnits = 'Hz'; % or 'rad/second', 'rpm', etc.
bode(G_dcm,options,{10,1e7})
grid on

%%

%Coefficient of the Compensator selected (PID transfer function)
fc=110e3; %crossover frequency
Ppm=60;  %phase margin
Pps=-100; %phase of the power stage at crosseover frequency fc
Pboost=-90+Ppm-Pps; %phase boost necessary to have a closed loop phase larger than th phase margin
Kboost=tand(45+(Pboost/2)); % k factor
fz=fc/Kboost;   % freqeuncy of the zero
wz=2*pi*fz;
fp=fc*Kboost;   %frequency of the pole
wp=2*pi*fp;
%%
%gains
Gpsd=-16.8;%gain of the power stage at corssover frequency in decibel
Gps=10^(abs(Gpsd/20));%gain of the power stage at corssover frequency
Gc=1/Gps;%gain of the compensator at corssover frequency
Kc=2*pi*fz/Gps; %gain of the compensator
%%
%%transfer function of the compensator
%%1 zero, 1 pair zero-pole
G1=Kc*tf([1/wz 1],[1 0]);
G2=tf( [0 1],[1/wp 1]);
Gc=G1*G2 % transfer function of the compensator
%bodw plot of teh compensator
options = bodeoptions;
options.FreqUnits = 'Hz'; % or 'rad/second', 'rpm', etc.
bode(Gc,options,{10,1e7})]