Code not running/plotting

% input parameters
Ta=15+273.15;
ha=5;
Tr=-7+273.15;
Tw=10+273.15;
Hi=10;
Bi=20;
Li=60;
Kins=0.045;
xbar=Li/Hi;
ybar=Bi/Hi;
er=0.05;
ep=0.95;
sigma=5.670e-8;
Tins=-5+273.15;
syms tins;

% define eqns
Aw1=Li*Bi;
Aw2=Hi*Bi;
Aw3=Li*Hi;
Awt=2*Aw3+2*Aw2
Rconv1=1/(ha*(2*Aw3+2*Aw2))
Rconv2=1/(ha*Aw1)
Rconv3=1/(ha*Aw1)
Rcond = tins/(Kins*Aw1)
F23=(2/(pi*xbar*ybar))*((log((((1+xbar^2)*(1+ybar^2))/(1+xbar^2+ybar^2))^0.5))+(xbar*(1+ybar^2)^(0.5))*atan((xbar)/((1+ybar^2)^0.5))+(ybar*(1+xbar^2)^0.5)*atan((ybar)/((1+xbar^2)^0.5))-xbar*atan(xbar)-ybar*atan(ybar));
Rrad2=1/(Aw1*F23)
F32=F23;
F31=1-F32;
Rrad3=1/(Aw1*F31)
F13=(Aw1*F31)/Awt
F12=1-F13
Rrad1=1/(Awt*F12)
Rrad4=(1-ep)/(ep*Aw3);


%CALCULATE NET q
syms Tc;
Eb3=sigma*Tc^4;
Eb1=sigma*Tw^4
Eb2=sigma*Tr^4;
Qrad=ep*Eb3*Aw1-ep*Eb1*Awt*F13-ep*Eb2*Aw1*F23;
Qconv=ha*Aw1*(Tc-Ta);
Qcond=(Tc-Tins)/(Rcond);


% Plotting T vs tins
Func = 0 ==  -Qconv -Qcond  -Qrad;
Tc = solve(Func,Tc);
fplot(Tc,[0.1,1])

How can I plot this code?

Answer 1

[Walter Roberson]

% input parameters
Ta=15+273.15;
ha=5;
Tr=-7+273.15;
Tw=10+273.15;
Hi=10;
Bi=20;
Li=60;
Kins=0.045;
xbar=Li/Hi;
ybar=Bi/Hi;
er=0.05;
ep=0.95;
sigma=5.670e-8;
Tins=-5+273.15;
syms tins;

% define eqns
Aw1=Li*Bi;
Aw2=Hi*Bi;
Aw3=Li*Hi;
Awt=2*Aw3+2*Aw2

Awt = 1600

Rconv1=1/(ha*(2*Aw3+2*Aw2))

Rconv1 = 1.2500e-04

Rconv2=1/(ha*Aw1)

Rconv2 = 1.6667e-04

Rconv3=1/(ha*Aw1)

Rconv3 = 1.6667e-04

Rcond = tins/(Kins*Aw1)

Rcond = 

F23=(2/(pi*xbar*ybar))*((log((((1+xbar^2)*(1+ybar^2))/(1+xbar^2+ybar^2))^0.5))+(xbar*(1+ybar^2)^(0.5))*atan((xbar)/((1+ybar^2)^0.5))+(ybar*(1+xbar^2)^0.5)*atan((ybar)/((1+xbar^2)^0.5))-xbar*atan(xbar)-ybar*atan(ybar));
Rrad2=1/(Aw1*F23)

Rrad2 = 0.0015

F32=F23;
F31=1-F32;
Rrad3=1/(Aw1*F31)

Rrad3 = 0.0018

F13=(Aw1*F31)/Awt

F13 = 0.3418

F12=1-F13

F12 = 0.6582

Rrad1=1/(Awt*F12)

Rrad1 = 9.4962e-04

Rrad4=(1-ep)/(ep*Aw3);


%CALCULATE NET q
syms Tc;
Eb3=sigma*Tc^4;
Eb1=sigma*Tw^4

Eb1 = 364.4595

Eb2=sigma*Tr^4;
Qrad=ep*Eb3*Aw1-ep*Eb1*Awt*F13-ep*Eb2*Aw1*F23;
Qconv=ha*Aw1*(Tc-Ta);
Qcond=(Tc-Tins)/(Rcond);


% Plotting T vs tins
Func = 0 ==  -Qconv -Qcond  -Qrad;
Tc = solve(Func,Tc, 'MaxDegree', 4);

tiledlayout('flow');
for K = 1 : 4;
    TcfunR = matlabFunction(real(Tc(K)));
    TcfunI = matlabFunction(imag(Tc(K)));
    
    nexttile();
    fplot({TcfunR, TcfunI}, [0.1, 1]); legend({'real', 'imag'}); title(string(K));
end

%fplot(Tc,[0.1,1])

[Valerie]

Thank you!