代码已出，最速下降法参数拟合，求改

tic
clc;clear;
syms x1 x2 x3 x4
t=[5 10 15 20 30 45 60 90 120];
p=[1;2;3;5;10;20;50;100];
i=[1.306 1.036 0.880 0.723 0.564 0.435 0.356 0.262 0.209
1.714 1.356 1.127 0.960 0.749 0.570 0.470 0.357 0.291
1.952 1.544 1.272 1.097 0.856 0.650 0.536 0.413 0.338
2.253 1.780 1.454 1.270 0.992 0.749 0.620 0.483 0.398
2.660 2.101 1.701 1.504 1.177 0.885 0.733 0.579 0.480
3.068 2.422 1.949 1.739 1.361 1.021 0.847 0.674 0.561
3.607 2.846 2.276 2.049 1.605 1.200 0.997 0.800 0.669
4.014 3.166 2.523 2.284 1.789 1.336 1.110 0.896 0.750];
for m=1:8
for n=1:9
F(m,n)=(x1.*(1+x2.*log(p(m,1)))/(t(1,n)+x3).^x4-i(m,n)).^2;
end
end
Q=sum(sum(F))
dQ1=diff(Q,x1);dQ2=diff(Q,x2);dQ3=diff(Q,x3);dQ4=diff(Q,x4);%对函数进行求一阶导
DQ=[dQ1;dQ2;dQ3;dQ4];
dQ11=diff(dQ1,x1);dQ12=diff(dQ1,x2);dQ13=diff(dQ1,x3);dQ14=diff(dQ1,x4);
dQ21=diff(dQ2,x1);dQ22=diff(dQ2,x2);dQ23=diff(dQ2,x3);dQ24=diff(dQ2,x4);
dQ31=diff(dQ3,x1);dQ32=diff(dQ3,x2);dQ33=diff(dQ3,x3);dQ34=diff(dQ3,x4);
dQ41=diff(dQ4,x1);dQ42=diff(dQ4,x2);dQ43=diff(dQ4,x3);dQ44=diff(dQ4,x4);
%这里进行求函数二阶导数
DQQ=[dQ11,dQ12,dQ13,dQ14;dQ21,dQ22,dQ23,dQ24;dQ31,dQ32,dQ33,dQ34;dQ41,dQ42,dQ43,dQ44];
%海赛矩阵
disp('x1 x2 x3 x4 分别代表 A c b n 取初值[20 0.9 10 1]')
x=input('请输入x1 x2 x3 x4的初始值为x=[x1 x2 x3 x4],x=:');
x=x';
%矩阵转置
E=input('请输入你所要求的最速下降法的精度数E=:');
%这里进行一些相关初始值的计算
T=[];
d=T;
TH=[];
%给d,t赋空矩阵
T(:,1)=subs(DQ,[x1,x2,x3,x4],[x(1),x(2),x(3),x(4)])
% T(:,1)指数组的第一列 res=subs(es,old,new)用new置换es中的old后产生res
TH=subs(DQQ,[x1,x2,x3,x4],[x(1),x(2),x(3),x(4)])
GG(1)=subs(Q,[x1,x2,x3,x4],[x(1),x(2),x(3),x(4)])
%这里进行一些相关初始值的计算
disp('由于你输入的初始值，这里将开始最速下降法搜寻：');
% 显示内容
for k=1:1000
d(:,1)=-T(:,1) %d(k)是x(k+1)=x(k)+R(k)*d(k)
%d(k)为计算搜索方向即负梯度方向
R(1)=(T(:,1)'*T(:,1))/(T(:,1)'*TH*T(:,1))
%R(k)为近似最佳步长
TH=subs(DQQ,[x1,x2,x3,x4],[x(1,k),x(2,k),x(3,k),x(4,k)])
T(:,k)=subs(DQ,[x1,x2,x3,x4],[x(1,k),x(2,k),x(3,k),x(4,k)])
d(:,k+1)=-T(:,k)
R(k)=(T(:,k)'*T(:,k))/(T(:,k)'*TH*T(:,k))
GG(k)=subs(Q,[x1,x2,x3,x4],[x(1,k),x(2,k),x(3,k),x(4,k)])
%while k>=2
%while GG(k)>=GG(k-1)
% R(k)=R(k)/2;
%end
% end
if norm(T(:,k))<E %norm（） 范-2误差
disp(' ');
disp('QX就是最速下降法的解 ')
QX=subs(Q,[x1,x2,x3,x4],[x(1,k),x(2,k),x(3,k),x(4,k)])
disp(' ');
disp('对应的x值为 ');
x(:,k)
break;
end
x(:,k+1)=x(:,k)+R(k)*d(:,k)
end
toc
各为兄弟姐妹，谁知道我编的这个哪里错了啊？
TH=subs(DQQ,[x1,x2,x3,x4],[x(1),x(2),x(3),x(4)]) 这句运行起来时间非常长。
有谁知道我这应该怎么改啊，希望好心人帮忙一下