Computes the inverse Laplace transform of a complex function.

The routine INLAP computes the inverse Laplace transform of a complex-valued function. Recall that if f is a function that vanishes on the negative real axis, then we can define the Laplace transform of f by

It is assumed that for some value of s the integrand is absolutely integrable.

The computation of the inverse Laplace transform is based on applying the epsilon algorithm to the complex Fourier series obtained as a discrete approximation to the inversion integral. The initial algorithm was proposed by K.S. Crump (1976) but was significantly improved by de Hoog et al. (1982). Given a complex-valued transform F(s) = L[f](s), the trapezoidal rule gives the approximation to the inverse transform

This is the real part of the sum of a complex power series in z = exp(iπ t/T), and the algorithm accelerates the convergence of the partial sums of this power series by using the epsilon algorithm to compute the corresponding diagonal Pade approximants. The algorithm attempts to choose the order of the Pade approximant to obtain the specified relative accuracy while not exceeding the maximum number of function evaluations allowed. The parameter α is an estimate for the maximum of the real parts of the singularities of F, and an incorrect choice of α may give false convergence. Even in cases where the correct value of α is unknown, the algorithm will attempt to estimate an acceptable value. Assuming satisfactory convergence, the discretization error E := g - f satisfies

It follows that if ∣f(t)∣ ≤ Meβt, then we can estimate the expression above to obtain (for 0 ≤ t ≤ 2T)

We invert the Laplace transform of the simple function (s - 1)-2 and print the computed answer, the true solution and the difference at five different points. The correct inverse transform is xex.