matlab ode解微分方程组，能否将微分方程组里的方程用符号表示？

代码如下：

[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])

function dcdzh2h2sco2_centri = dcdzh2h2sco2_centri(z,c)
dcdzh2h2sco2_centri = zeros(3,1);
dcdz_1 = (27*((c(3)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(2)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - 1))/(1000*(10*z + 293)) + 1/((4157*z)/1700000 + 1218001/17000000);
dcdzh2h2sco2_centri(1) = dcdz_1;
% dcdzh2h2sco2_centri(1) = (27*((c(3)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(2)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - 1))/(1000*(10*z + 293)) + 1/((4157*z)/1700000 + 1218001/17000000);
dcdzh2h2sco2_centri(2) = 1/((4157*z)/2200000 + 1218001/22000000) - (27*((c(1)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - (c(3)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + 1))/(1000*(10*z + 293));
dcdzh2h2sco2_centri(3) = 1/((4157*z)/800000 + 1218001/8000000) - (27*((c(2)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(1)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + 1))/(1000*(10*z + 293));
end

以上代码是可以运行的，在微分方程组里，先把第一个微分方程dcdzh2h2sco2_centri(1)右面的定义为dcdz_1，然后再让dcdzh2h2sco2_centri(1)等于dcdz_1，和直接定义dcdzh2h2sco2_centri(1)（被注释掉的那部分）结果一样；

可是微分方程是由前面符号推导得出的，如果直接运行下面的代码：

[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])

function dcdzh2h2sco2_centri = dcdzh2h2sco2_centri(z,c)

R = 8.314;
g = 1000;
Mi_h2s = 34;
Mi_co2 = 44;
Mi_ch4 = 16;
syms T z c_ch4 c_co2 c_h2s
%%
% entrance conditions:T = 293K;p = 0.10325Pa;
% entrance mass fraction:h2s:co2 = 0.156:0.078
% c_ch4(0)=37.063--3
% c_co2(0)=1.3932--2
% c_h2s(0)=3.6029--1
Tz = 293 + 10*z;%温度随深度的变化
alpha_1_2 = log((3.6029/(1.3932))*((c_co2)/c_h2s))/(log(Tz/293.15));
alpha_1_3 = log((3.6029/(37.063))*((c_ch4)/c_h2s))/(log(Tz/293.15));
alpha_2_3 = log((1.3932/(37.063))*((c_ch4)/c_co2))/(log(Tz/293.15));
alpha_h2sT = (c_co2/(c_ch4+c_co2+c_h2s))*(-alpha_1_2)+((c_ch4)/(c_ch4+c_co2+c_h2s))*(-alpha_1_3);
alpha_co2T = (c_ch4/(c_ch4+c_co2+c_h2s))*(-alpha_2_3)+((c_h2s)/(c_ch4+c_co2+c_h2s))*(alpha_1_2);
alpha_ch4T = (c_h2s/(c_ch4+c_co2+c_h2s))*(alpha_1_3)+((c_co2)/(c_ch4+c_co2+c_h2s))*(alpha_2_3);
alpha_h2s = subs(alpha_h2sT,T,Tz);
alpha_co2 = subs(alpha_co2T,T,Tz);
alpha_ch4 = subs(alpha_ch4T,T,Tz);
Hi_h2s = 1/((Mi_h2s*g)/(R*Tz));
Hi_co2 = 1/((Mi_co2*g)/(R*Tz));
Hi_ch4 = 1/((Mi_ch4*g)/(R*Tz));
dcdz_h2s_1 = -(0.027*((1 + alpha_h2s)/Tz)-(1/Hi_h2s))
dcdz_co2_1 = (-(0.027*((1 + alpha_co2)/Tz)-(1/Hi_co2)))
dcdz_ch4_1 = (-(0.027*((1 + alpha_ch4)/Tz)-(1/Hi_ch4)))
dc_11 = char(dcdz_h2s_1)
dc_12 = char(dcdz_co2_1)
dc_13 = char(dcdz_ch4_1)
dc_21 = strrep(dc_11,'c_h2s','c(1)');
dc_21 = strrep(dc_21,'c_co2','c(2)');
dc_21 = strrep(dc_21,'c_ch4','c(3)')
dc_22 = strrep(dc_12,'c_h2s','c(1)');
dc_22 = strrep(dc_22,'c_co2','c(2)');
dc_22 = strrep(dc_22,'c_ch4','c(3)')
dc_23 = strrep(dc_13,'c_h2s','c(1)');
dc_23 = strrep(dc_23,'c_co2','c(2)');
dc_23 = strrep(dc_23,'c_ch4','c(3)')
% global dc_31
dc_31 = str2sym(dc_21)
%%
dcdzh2h2sco2_centri = zeros(3,1);
dcdz_1 = (27*((c(3)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(2)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - 1))/(1000*(10*z + 293)) + 1/((4157*z)/1700000 + 1218001/17000000);
dcdzh2h2sco2_centri(1) = dcdz_1;
% dcdzh2h2sco2_centri(1) = (27*((c(3)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(2)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - 1))/(1000*(10*z + 293)) + 1/((4157*z)/1700000 + 1218001/17000000);
dcdzh2h2sco2_centri(2) = 1/((4157*z)/2200000 + 1218001/22000000) - (27*((c(1)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) - (c(3)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + 1))/(1000*(10*z + 293));
dcdzh2h2sco2_centri(3) = 1/((4157*z)/800000 + 1218001/8000000) - (27*((c(2)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + (c(1)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(3) + c(2) + c(1))) + 1))/(1000*(10*z + 293));
end

这里的dcdzh2h2sco2_centri(1)就是dc_21，会提示：
从 sym 转换为 double 时出现以下错误:
Unable to convert expression into double array.

出错 dcdzh2h2sco2_centri (line 51)
dcdzh2h2sco2_centri(1) = dc_31;

出错 odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

出错 ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 a_0320_centrifuge_100g_293k (line 54)
[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])
那么该如何在微分方程组里把某个方程的右端用一个符号表达式代替，这个符号表达式该如何转化为double来进行被识别？

你好，建议这么写，绝对很舒服（有帮助还望题主给个宝贵的采纳支持一下呢，谢谢啦）

[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])
 
function dcdzh2h2sco2_centri = dcdzh2h2sco2_centri(z,c)
 
R = 8.314;
g = 1000;
Mi_h2s = 34;
Mi_co2 = 44;
Mi_ch4 = 16;
c = num2cell(c(:)');                    % 关键在这两行
[c_h2s, c_co2, c_ch4] = c{:}; % 关键在这两行
syms T
%%
% entrance conditions:T = 293K;p = 0.10325Pa;
% entrance mass fraction:h2s:co2 = 0.156:0.078
% c_ch4(0)=37.063--3
% c_co2(0)=1.3932--2
% c_h2s(0)=3.6029--1
Tz = 293 + 10*z;%温度随深度的变化
alpha_1_2 = log((3.6029/(1.3932))*((c_co2)/c_h2s))/(log(Tz/293.15));
alpha_1_3 = log((3.6029/(37.063))*((c_ch4)/c_h2s))/(log(Tz/293.15));
alpha_2_3 = log((1.3932/(37.063))*((c_ch4)/c_co2))/(log(Tz/293.15));
alpha_h2sT = (c_co2/(c_ch4+c_co2+c_h2s))*(-alpha_1_2)+((c_ch4)/(c_ch4+c_co2+c_h2s))*(-alpha_1_3);
alpha_co2T = (c_ch4/(c_ch4+c_co2+c_h2s))*(-alpha_2_3)+((c_h2s)/(c_ch4+c_co2+c_h2s))*(alpha_1_2);
alpha_ch4T = (c_h2s/(c_ch4+c_co2+c_h2s))*(alpha_1_3)+((c_co2)/(c_ch4+c_co2+c_h2s))*(alpha_2_3);
alpha_h2s = subs(alpha_h2sT,T,Tz);
alpha_co2 = subs(alpha_co2T,T,Tz);
alpha_ch4 = subs(alpha_ch4T,T,Tz);
Hi_h2s = 1/((Mi_h2s*g)/(R*Tz));
Hi_co2 = 1/((Mi_co2*g)/(R*Tz));
Hi_ch4 = 1/((Mi_ch4*g)/(R*Tz));
dcdz_h2s_1 = -(0.027*((1 + alpha_h2s)/Tz)-(1/Hi_h2s));
dcdz_co2_1 = (-(0.027*((1 + alpha_co2)/Tz)-(1/Hi_co2)));
dcdz_ch4_1 = (-(0.027*((1 + alpha_ch4)/Tz)-(1/Hi_ch4)));

dcdzh2h2sco2_centri = zeros(3,1);
dcdzh2h2sco2_centri(1) = dcdz_h2s_1;
dcdzh2h2sco2_centri(2) = dcdz_co2_1;
dcdzh2h2sco2_centri(3) = dcdz_ch4_1;
end