Greetings. From the link above, the description of Post's formula states,

"As can be seen from the formula, the need to evaluate derivatives of arbitrarily high orders renders this formula impractical for most purposes."

With all of that said, using the Bromwich Integral might be a better method for higher orders of "s".

> # Example 1 with small value of exponent(3):
> restart:
> L(s):=2/s^3;
> ls:=exp(s*t)*L(s);
> denominator:=factor(denom(ls));
> p:=solve(denom(ls)=0,s);
> res:=px->residue(ls,s=p[px]);
> res1:=res(1);
> evalf(subs(t=10.0,res1));

The answer matches your procedure above.

> # Example 2 with larger value of exponent(5):
> restart:
> L(s):=2/s^5;
> ls:=exp(s*t)*L(s);
> denominator:=factor(denom(ls));
> p:=solve(denom(ls)=0,s);
> res:=px->residue(ls,s=p[px]);
> res1:=res(1);
> evalf(subs(t=10.0,res1));

Here we get an answer that matches what intrans(invlaplace) yields with the substitution of t=10. Because of the higher order of s, the procedure you have listed above does not return a numerical value.

Regards,
Georgios Kokovidis

P.S. Welcome back to Maple. I missed your posts. Most of what I know in Maple I learned from studying examples like yours. Thanks.

num_invLaplace_1:=proc(F,tau)
# Stehfest algorithm, according to Kou
#option remember;
local oldDigits, w,g, result, n,nMax, r,k, t;

nMax:=20;
oldDigits:=Digits: Digits:=oldDigits + max(80, 2*nMax);

w:= unapply( (-1)^(n-k)*k^n/k!/(n-k)!, k,n);
g:= unapply( ln(2)/t*(2*k)!/k!/(k-1)!*Sum( (-1)^r*binomial(k,r)*F((k+r)*ln(2)/t)  ,r=0..k), k,t):

t:=evalf(tau);
result:= 0;

k:=1;
while ( k LE nMax ) do
  result:= evalf(result + w(k, nMax)*g(k,t));
  k := k+1;
end do;

return evalf[oldDigits](result);
end proc;

Read 'LE' as less or equal (i never understand how to input that here ...)

some test functions:

#F:= s -> 1 / (s + 0.5);  F:= s -> 1/((s^2-1)^(1/2)); F:= s -> s^2/(s^2+2^2)^(3/2);
F:= s -> 2/s^8;
inttrans[invlaplace](F(s), s, T):
f:=unapply(%,T); #plot( f(t), t=0..8);

                                       2
                            F := s -> ----
                                        8
                                       s

                                           7
                         f := T -> 1/2520 T
some test values:

tTst:=10;
num_invLaplace_1(F,tTst); #num_invLaplace(F,tTst); 
evalf(f(tTst)); 
`approximation error` = evalf(%-%%);
                              tTst := 10

                           3968.2539682541
                           3968.2539682540
                                                -9
                   approximation error = -0.1 10

The other test functions involve Bessel functions and need more steps.

One uses an argument of Abate & Whitt to transform the inversion problem
to the form Int( g(s)*cos(t*s) , s = 0 .. infinity), where g is real and
analytic on the positive Reals and some real constant, see the file below,
which is easily chosen depending on the function to be inverted (it does
not depend on the point in which one wants to evaluate the Laplace inverse).

These integrals with oscillating factor can be computed efficiently according
to Mori & Ooura (the latter provides C and Fortran code, translate that to
Maple and mildly increase precision for the nasty inversion problem).

Download 102_A_Way_to_Compute_the_Numerical_Laplace_Inversion.zip
View file details

Have you had a look at Zakian's Method of numericaly inverting Laplace transforms, quick, simple but does have its problems

Please replace this text with the link to your file.
The link can be found in the File Manager

If there is not a problem with the uploading mechcanism, perhaps someone can write a procedure on how to do it, as I have tried to upload this worksheet any number of times!