请问：如何利用matlab实现Gaver-stehfest 算法的拉普拉斯变换？

我写的程序，可是效果很差，请高手指正

format long
n=9;
t=1;
sum1=0;
p=0;
for k=1:n
    d1=(-1)^(n-k)*k^n/(factorial(k)*factorial(n-k))*log(2)/t...
        *factorial(2*k+4)/(factorial(k+2)*factorial(k+1));
    for m=0:k+2
        f=@(x)x*0.02+0.5*x^2*0.05^2+2*(0.5*30/(30-x)+0.5*30/(30+x)-1)-(m+k+2)*log(2)/t;
        x1=fzero(f,-60);
        x2=fzero(f,-20);
        a1=(30+x2)/30*(-x1)/(x2-x1);
        a2=(-x1-30)/30*(-x2)/(x2-x1);
        d3=t/((m+k+2)*log(2))*(a1*exp(0.9*(-x2))+a2*exp(0.9*(-x1)));
        d2=(-1)^m*factorial(k+2)/(factorial(m)*factorial(k+2-m));
        sum1=sum1+d2*d3;
    end
    p=p+d1*sum1;
end
p

[此贴子已经被作者于2008-7-5 22:34:54编辑过]