I have a mathematica code for solving numerical inverse Laplace transform (Credit to Mr. Patrick O. Kano), but sadly the code is using `Nintegrate`, and i want to use Simpson's rule with adjustable stepsize, and function. By the way my numerical ILT is based on Talbot's method. For more details about Talbot's method, parametrization and its contour, i've mentioned it on my previous question here:
[Need help to plot Talbot's Contour (Newbie question)][1]

Hope this doesn't sound complicated. So here is my Mathematica code for numerical ILT.

    Timeval = 2;
    Rval = 2;
    Flap[s_] = 1/(s - 1);
    Tfunexact[t_] = Exp[1*t];
    Valexact = N[Tfunexact[Timeval], 10]
    
    STalbot[r_, a_] = r*a*Cot[a] + I*r*a;
    dsda[r_, a_] = (1 + I*(a + Cot[a]*(a*Cot[a] - 1)));
    TimeDfun[r_, t_] = 
      Rval/(Pi)*
       NIntegrate[
        Exp[STalbot[r, a]*t]*Flap[STalbot[r, a]]*dsda[r, a], {a, 0, Pi}, 
        WorkingPrecision -> 10];
    
    {Timeval, Approxval} = Timing[TimeDfun[Rval, Timeval]]
    RelError = Abs[Approxval - Valexact]/Valexact

Notice the `Nintegrate` there, i want to replace it with Simpson's rule, and here is my attempt:

    Clear[a, b, n]
    simpson[r_, 
      t_] := (Rval/(3 n)) ( 
       r*Exp[r t]*Flap[r] +  
        2 Sum[Exp[STalbot[r, 2 i \[Theta]]*t]*Flap[STalbot[r, \[Theta]]]*
           dsda[r, \[Theta]]] +  
        4 Sum[Exp[STalbot[r, (2 i - 1) \[Theta]]*t]*
           Flap[STalbot[r, (2 i - 1) \[Theta]]]*
           dsda[r, (2 i - 1) \[Theta]]]   )
        
    f[s_] := 1/(s - 1)
    TableForm[
     Table[{n, N[simpson[0, \[Pi], n]]}, {n, 10, 100, 10}, 
      TableHeadings -> {{}, {"n", "s_n"}}]]

Can you please help me to accomplish this. I was failed to perform the last code. Please comment if there's something unclear about my post. Thanks in advance!


  [1]: https://mathematica.stackexchange.com/questions/240487/need-help-to-plot-talbots-contour-newbie-question