急求96年数学建模题A题的解答

%保存为M文件，文件名为"smax1.m"，在主窗口中输入"[ss3,ss4,st]=smax1()"命令来运行！

%k为捕捞强度，通过改变k来改变结果。

％本程序还有很多不足之处，请大家多多指教！！

function [ss3,ss4,st]=smax1()
s10=122e9;
s20=29.7e9;
s30=10.1e9;
s40=3.29e9;
s0=[s10,s20,s30,s40];
a1=0.8;
ss3=0;
ss4=0;
a=1.109e5;
t1=1;
syms t;
d=10;
k=10;
h=1e5;
while d>=0
sb1=s10*exp(-a1*(2/3)*t1);
sb2=s20*exp(-a1*(2/3)*t1);
sb3=s30*exp(-(a1+0.42*k)*(2/3)*t);
sb4=s40*exp(-(a1+k)*(2/3)*t);

%sb3q=s30*exp(-(a1+0.42*k)*(2/3)*t1);

s1=sb1*exp(-(1/3)*t1);
s2=sb2*exp(-(1/3)*t1);
%s3=sb3*exp(-(1/3)*t);
%s4=sb4*exp(-(1/3)*t);

s3q=s30*exp((-(a1+0.42*k)*(2/3)-(1/3))*t1);
s3=s30*exp((-(a1+0.42*k)*(2/3)-(1/3))*t);
s4=s40*exp((-(a1+k)*(2/3)-(1/3))*t);

sp3=quad8('10.1e9*exp((-(0.8+0.42*10)*(2/3)-(1/3))*t)/4',0,1/3)/(1/3);
sp4=quad8('3.29e9*exp((-(0.5+10)*(2/3)-(1/3))*t)/4',0,1/3)/(1/3);
n=sp4*a+sp3*a*0.5;
s11=n*1.22e11/(1.22e11+n);
s21=s1;
s31=s2;
s41=s3q;

c=[s11,s21,s31,s41];

if abs(s11-s10) ss3=s31;
ss4=s41;
st=17.86*s30*exp(-(0.42*k)*(2/3)*t1)+22.99*s40*exp(-k*(2/3)*t1);
d=-1;
else
s10=s11;
s20=s21;
s30=s31;
s40=s41;
d=d-1;
end
end

中科大数学建模网站上有1993-2002年历年的解答，自己去下