我审查完公式，又从新编辑，希望你能再提点一下

我的想法是 解四个方程，a，b是参数，t(1)_t(4)为变量，t(1)_t(4)分别随a，b变化，生成三维图，

t0=[0;7;0];

a=0:5;

b=0:5;

[A,B]=meshgrid(a,b)

p1=zeros(6,6);

p2=zeros(6,6);

p3=zeros(6,6);

p4=zeros(6,6);

for i=1:length(a)

for j=1:length(b)

p(:,i)=fsolve(@(t)SH_2(t,a(i),b(j)),t0);

p1(i,j)=real(sol(1));

p2(i,j)=real(sol(2));

p3(i,j)=real(sol(3));

p4(i,j)=real(sol(4));

end

end

figure

subplot(2,2,1)

surf(A,B,p1); shading interp; colorbar

subplot(2,2,2)

surf(A,B,p2); shading interp; colorbar

subplot(2,2,3)

surf(A,B,p3); shading interp; colorbar

subplot(2,2,4)

surf(A,B,p4); shading interp; colorbar

function M = SH_2( t,a,b )

y(1)=exp(20.37788-4.60517*t(2))+exp(30.14084+log(t(1))-4.60517*t(2))+exp(38.40712+2*log(t(1))-4.60517*t(2))+exp(45.2458+3*log(t(1))-4.60517*t(2))+exp(50.31148+4*log(t(1))-4.60517*t(2))+exp(48.99901+5*log(t(1))-4.60517*t(2))+exp(2.32561-2.30259*t(2))+exp(-15.45035)+exp(-44.18661+2.30259*t(2))+exp(-74.30442+4.60517*t(2))+exp(15.24311-4.60517*t(2))+exp(15.24311-4.60517*t(2))+exp(22.03574+3*log(t(1))-2.30259*t(2))+exp(-8.77285+2*log(t(1)))+exp(-40.04195+log(t(1))+2.30259*t(2))-t(4);

y(2)=exp(log(t(1)))+exp(22.56533+log(t(1))-2.30259*t(2))+exp(30.14084+log(t(1))-4.60517*t(2))+2*exp(38.40712+2*log(t(1))-4.60517*t(2))+3*exp(45.2458+3*log(t(1))-4.60517*t(2))+4*exp(50.31148+4*log(t(1))-4.60517*t(2))+5*exp(48.99901+5*log(t(1))-4.60517*t(2))+3*exp(22.03574+3*log(t(1))-2.30259*t(2))+2*exp(-8.77285+2*log(t(1)))+exp(-40.04195+log(t(1))+2.30259*t(2))-a-2*b;

y(3)=exp(22.56533+log(t(1))-2.30259*t(2))+exp(-2.30259*t(2))-exp(-32.65066+2.30259*t(2))+2*exp(20.37788-4.60517*t(2))-2*exp(log(b))+2*exp(30.14084+log(t(1))-4.60517*t(2))+2*exp(38.40712+2*log(t(1))-4.60517*t(2))+2*exp(45.2458+3*log(t(1))-4.60517*t(2))+2*exp(50.31148+4*log(t(1))-4.60517*t(2))+2*exp(48.99901+5*log(t(1))-4.60517*t(2))+exp(2.32561-2.30259*t(2))-exp(-44.18661+2.30259*t(2))-2*exp(-74.30442+4.60517*t(2))+2*exp(15.24311-4.60517*t(2))+exp(22.03574+3*log(t(1))-2.30259*t(2))-exp(-40.04195+log(t(1))+2.30259*t(2));

M=[y(1);y(2);y(3)];

end