求高手帮忙！！这个程序运行时间能缩短吗？

clear
Pij=[1.0705 0.6828 0.6345 0.4370 0.4967 0.4450] %各段的平均载荷矩阵
h=[9.84 19.68 29.52 39.36 44.34 49.34] %各段的高度矩阵
for i=1:6
for j=1:6
Spij(i,j)=exp(-abs(h(i)-h(j))/60).*Pij(i).*Pij(j)
end
end
goal=chol(Spij)
goal=goal'
N=4001
afa=2.205
Wf=0.109
kexi=1.485
dw=0.01256
V=26.8
x=sqrt(-1)
fai=2*pi*rand(1,N)
for time=1:4001
t=(time-1)*0.05
for j=1:6
sum1=0
for m=1:j
sum2=0
for k=1:N
w=(k-1)*dw
F1=sqrt(afa*pi*V/900)*w*x
F2=(Wf^2-w^2)+2*kexi*Wf*w*x
GS=F1/F2
GD(j,m)=sqrt(2*dw)*goal(j,m)*abs(GS)
if k<=1
theta=0.0
end
if k>1
theta=atan2(imag(GS),real(GS))
end
Result=GD(j,m)*cos(w*t+theta+fai(k))
sum2=sum2+Result
end
xie(j,m)=sum2
sum1=sum1+xie(j,m)
end
ft(time,j)=sum1
end
ok(time)=t
end
for l=1:4001
tu1(l)=ft(l,1)
tu2(l)=ft(l,2)
tu3(l)=ft(l,3)
tu4(l)=ft(l,4)
tu5(l)=ft(l,5)
tu6(l)=ft(l,6)
end
运行了半个小时还在进行，然后忽然自己退出了啊，悲剧的，求高手指点，这个程序有没有缩短运行时间的办法？另外我用41替换过4001，是能出结果的。程序应该不是死循环吧。。。