syms t
y6=107163*(0.24*(51749*t)^0.2+0.35)^2
y5=0.0526*log(0.0001955*t*y6+1)
y6=matlabFunction(y6);
y5=matlabFunction(y5);
odefunc=@(t,y)[111.732*y5(t)*y(2)-(111.732*y5(t)*y(1))/47500;
(111.732*y5(t)*y(1)+0.0000936*y(3))-(y(2)*(111.732*y5(t)+0.00243))/9175;
(0.00243*y(2)+0.326*111.732*y5(t))*y(4)-(y(3)*(0.0000936+0.326*111.732*y5(t))/(3.9765*10^12 ));
(1+0.326*111.732*y5(t)*y(3))-(y(4)*0.326*111.732*y5(t))/1329500
]
tspan=[2,2.9]; %时间区间
t0=[ 0.0000000155705,0,0.000032,0];
[t,y]=ode45(odefunc,tspan,t0);
plot(t,y(:,1))
方程如果都能够显示表达的话，matlab求解也是挺方便的。mathematica之所以一直busy，应该跟时间区间的选取有关，时间区间大于2.9的时候，matlab也会给出警告
