I have used GWR for many years, with very good results. This is the illusttration how to use it:

(*http://library.wolfram.com/infocenter/MathSource/4738/*)GWR[F_,t_,M_: 32,precin_: 0]:=Module[{M1,G0,Gm,Gp,best,expr,τ=Log[2]/t,Fi,broken,prec},If[precin<=0,prec=21 M/10,prec=precin];
If[prec<=$MachinePrecision,prec=$MachinePrecision];
broken=False;
If[Precision[τ]<prec,τ=SetPrecision[τ,prec]];
Do[Fi[i]=N[F[i τ],prec],{i,1,2 M}];
M1=M;
Do[G0[n-1]=τ (2 n)!/(n! (n-1)!) Sum[Binomial[n,i] (-1)^i Fi[n+i],{i,0,n}];
If[Not[NumberQ[G0[n-1]]],M1=n-1;G0[n-1]=.;Break[]];,{n,1,M}];
Do[Gm[n]=0,{n,0,M1}];
best=G0[M1-1];
Do[Do[expr=G0[n+1]-G0[n];
If[Or[Not[NumberQ[expr]],expr==0],broken=True;Break[]];
expr=Gm[n+1]+(k+1)/expr;
Gp[n]=expr;
If[OddQ[k],If[n==M1-2-k,best=expr]];,{n,M1-2-k,0,-1}];
If[broken,Break[]];
Do[Gm[n]=G0[n];G0[n]=Gp[n],{n,0,M1-k}];,{k,0,M1-2}];
best]
SetAttributes[GWR,{Listable,NHoldAll}]

F[s_]:=12 (1-Exp[168 s]+168 Exp[180 s] s)/(-1+Exp[12 s]+Exp[168 s]-Exp[180 s]+2016 Exp[204 s] s^2);

t={10^-1,1,10,100,1000,10000,100000};

N[GWR[F,t,60],20]
N[GWR[F,t,80],20]
N[GWR[F,t,100],20]
N[GWR[F,t,120],20]
N[GWR[F,t,140],20]

You should increase the number of digits (the third calling parameter of GWR) until you get a satisfactory accuracy of the results. I notice that your original formula has one extra parenthesis, so that I am not sure if after removing it the formula for F is correct. But if it is, then I see that your transform may require a rather high value of the third parameter. Of course, instead of using a list of time values (t) you can just put a single value. Leslaw

Which method is the best (you already used Durbin)?

This is a difficult question and I have no competence and I am not an expert in this field.In my humble opinion, it depends on the given function what method is the best.You just have to experiment yourself.

See attached file.

Attachments:

InverseLaplace.nb

See attached file.

Attachments:

InverseLaplace ver2.nb

See attached file.

Attachments:

InverseLaplace ver3.nb

I had copied and pasted from your second line. Now it's better:

f[t_] = InverseLaplaceTransform[(12  E^(24 s) 
      (1 - E^(168   s) + 168   E^(180   s)   s))/
    (-1 + E^(12   s) + E^(168   s) -
      E^(180   s) + 2016   E^(204   s)   s^2), s, t];
{f[1.], f[10.], f[1000.], 16/19.}

{1., 0.999996, 0.842102, 0.842105}

I am no expert on Laplace transforms, I just tried out the built-in commands of Mathematica:

f[t_] = InverseLaplaceTransform[(12  (1 - E^(168  s) + 
       168  E^(180  s)  s))/(-1 + E^(12  s) + E^(168  s) - 
     E^(180  s) + 2016  E^(204  s)  s^2), s, t]

gives no closed-form output. I tried numerical values:

In[1986]:= {f[1.], f[10.], f[1000.], 16/19.}

Out[1986]= {1.27295*10^-111, -1.82483*10^-10, 0.8421, 0.842105}

which partially agree with what you know.

If I compute first term of asymptotic I have: $\theta (t-24)$

and looks like for: 0<t<24 is zero and for t > 24 is 1.

 func = (12 (1 - E^(168 s) + 168 E^(180 s) s))/(-1 + E^(12 s) + E^(
     168 s) - E^(180 s) + 2016 E^(204 s) s^2);
 g[t_] := InverseLaplaceTransform[func, s, N@t, Method -> "Durbin"];
 f = InverseLaplaceTransform[Asymptotic[func, {s, Infinity, 1}], s, t]
  Show[{Plot[f, {t, 0, 50}, PlotStyle -> Red], 
    ListLinePlot[Table[{t, g[t]}, {t, 1/10, 100, 1/3}]] // Quiet}, 
   PlotRange -> All]

Plot looks correct to me.

Regards M.I