matlab solve 带入数值,solve求解方程组怎么不自动带入系数的数值？

clear all

clc

G_w=6;

p_g=2500;

c_g=840;

H_g=0.0032;

h_a=2.8;

xi_g=0.1;

k_a=2.76;

H_a=0.0001;

p_pv=2330;

c_pv=916;

H_pv=0.002;

K_pv=150;

xi_TPT=0.27;

xi_pv=0.81;

r=0.95;

A_pv=0.048;

r1=0.038;

n_pv=0.15;

B_pv=0.0045;

T_n=25;

K_ei=2;

H_ei=0.001;

R_ei=H_ei/K_ei;

p_Al=2700;

c_Al=880;

K_Al=237;

H_Al=0.0002;

r2=0.315;

h_Al_a=11;

M_eva=0.137;

M_con=0.275;

c_pipe=384;

K_pipe=398;

D_con=0.015;

R_evacon=0.000041;

A_hp_pv=0.003*1.26;

u=G_w/(60000*9.6*10^(-4));

Re=0.015*u/(0.805*10^(-6));

Nu=0.027*Re^0.63*5.42^0.36*(5.42/4.31)^0.25;

K_w=0.618;

h_w_con=Nu*K_w/D_con;

A_hp_w=0.035*pi*D_con;

M_w=(G_w/60)/1000*1000;

m_w=0.072;

c_w=4183;

R_w_a=11;

M_tank=120;

T0=30;

b=5.6*10^(-8);

simtime=2*60;

dt=60;

G=800;

T_a=zeros(simtime/dt+1,1);            %室外空气温度定义为矩阵

T_sky=zeros(simtime/dt+1,1);          %有效天空温度

G=zeros(simtime/dt+1,1);              %太阳辐射强度

T_g=zeros(simtime/dt+1,1);            %玻璃板的温度

T_pv=zeros(simtime/dt+1,1);           %pv板的温度

T_Al=zeros(simtime/dt+1,1);           %铝片的温度

T_hp_con=zeros(simtime/dt+1,1);       %热管冷凝段温度

T_hp_eva=zeros(simtime/dt+1,1);       %热管蒸发段温度

T_w_box=zeros(simtime/dt+1,1);        %联箱中水的温度

T_w_tank=zeros(simtime/dt+1,1);       %水箱中水的温度

T_g(1,1)=T0;

T_pv(1,1)=T0;

T_Al(1,1)=T0;

T_hp_con(1,1)=T0;

T_hp_eva(1,1)=T0;

T_w_box(1,1)=T0;

T_w_tank(1,1)=T0;

for k=2:(simtime/dt+1)

syms Ta Tb Tc Td Te Tf Tg

eq1=(p_g*c_g*H_g)*(Ta-30)/dt-(h_a*(30-Ta)+xi_g*b*((9.07+273.15)^2+(30+273.15)^2)*(30+30+2*273.15)*(9.07-Ta)+(0.111*b*((30+273.15)^2+(30+273.15)^2)*(30+30+2*273.15)+k_a/H_a)*(Tb-Ta)+G*xi_g);

eq2=(p_pv*c_pv*H_pv)*(Tb-30)/dt-(800*xi_TPT*(1-r)+G*xi_pv*r-G*xi_pv*r*n_pv*(1-B_pv*(Tb-25))+(0.111*b*((30+273.15)^2-(30+273.15)^2)*(30+30+2*273.15)+k_a/H_a)*(Ta-Tb)+r1/R_ei*(Td-Tb)+(1-r1)/R_ei*(Tc-Tb));

eq3=(p_Al*c_Al*H_Al)*(Tc-30)/dt-((1-r1)/R_ei*(Tb-Tc)+(0.315-r1)/R_ei*(Td-Tc)+(30-Tc)*h_Al_a);

eq4=(M_eva*c_pipe)*(Td-30)/dt-((Te-Td)/R_evacon+(Tb-Td)*A_hp_pv/R_ei+(Tc-Td)*(0.008*pi*1.26-A_hp_pv)/R_ei);

eq5=(M_con*c_pipe)*(Te-30)/dt-((Td-Te)/R_evacon+12*A_hp_w*h_w_con*(Tf-Te));

eq6=M_w*c_w*(Tf-20)-(A_hp_w*12*h_w_con*(Te-Tf)+0.9*0.1*4*h_Al_a*(20-Tf));

eq7=M_tank*c_w*(Tg-20)/dt-(M_w*c_w*(Tf-20));

s=solve(eq1,eq2,eq3,eq4,eq5,eq6,eq7);

s=double([s.Ta s.Tb s.Tc s.Td s.Te s.Tf s.Tg]);

T_g(k,1)=s(1);

T_pv(k,1)=s(2);

T_Al(k,1)=s(3);

T_hp_eva(k,1)=s(4);

T_hp_con(k,1)=s(5);

T_w_box(k,1)=s(6);

T_w_tank(k,1)=s(7);

end