How I can plot this equation?

TB =3.8470608473110835049785363403119;
TI =15.448004476144693589582536852599;
n20av=1.1;
%confinement time [s]
tauE=1.2;
%internal radius [m]
a=2.0;
%plasma ellipticity
k=1.8;
%external radius [m]
R0=6.2;
%inverse aspect rdius
epsy=a/R0;
%plasma volume [m3]
Vp=2*pi*R0*pi*k*a^2;
%plasma power [W/m3]
%magnetic current [MA]
Imax=7.5;
%homic power [W/m3]

%bremmstraulung power [W/m3]
Sb1=6.14*(10^3)*(n20av^2)*TB.^(1/2);
%conduction power
Sk2=5.34*(10^4)*n20av*TI*1/tauE;

T1=0.1;
sigmav1=(10^(-6))*exp(-21.38*(T1.^(-0.2935))-25.20-7.101*(10^(-2))*T1+1.938*(10^(-4))*T1.^2+4.925*(10^(-6))*T1.^(3)-3.984*(10^(-8))*T1.^4);
sigmavn1=sigmav1/10^(-22);
Sa1=(2.31*10^(5))*(n20av^(2))*sigmavn1;
Sohm1=(1/Vp)*(5.6*10^(-2)*((1-1.31*epsy^(1/2)+0.46*epsy)).^-1)*(R0*Imax^(2))*((a^(2)*k*((T1).^(3/2))).^-1);
dwdt = Sa1 + Sohm1-Sb1-Sk2;

syms Tk
W_t=(5.34*10^(4))*n20av*Tk;
bil = W_t +dwdt ==0;
w=vpasolve(bil,Tk);

T=linspace(0,0.5,200);
sigmav1=(10^(-6))*exp(-21.38*(T.^(-0.2935))-25.20-7.101*(10^(-2))*T+1.938*(10^(-4))*T.^2+4.925*(10^(-6))*T.^(3)-3.984*(10^(-8))*T.^4);
sigmavn1=sigmav1/10^(-22);
Sa1=(2.31*10^(5))*(n20av^(2))*sigmavn1;
Sohm1=(1/Vp)*(5.6*10^(-2)*((1-1.31*epsy^(1/2)+0.46*epsy)).^-1)*(R0*Imax^(2))*((a^(2)*k*((T).^(3/2))).^-1);
dwdt = Sa1 + Sohm1-Sb1-Sk2;
bill = w+dwdt;

for T=linspace(1,16,400)
  
    if (TB< T) && T <= TI
         Saux = (1-((T-TB)/(TI-TB)));
         bill1= bill+Saux;
       
    elseif T>TI
        Saux=0;
    end
end

plot(T,bill1)

Answer 1

[Mark Sherstan]

plot(bill1)

You are currently getting an error as T is only 1x1 and bill1 is 1x200

[Franco Belletti]

thanks!