
function f = direct(F, M, t, at)

%  The function f = direct(F, M, t, at) inverts the Laplace
%  transform by a method based on a trapezoidal rule approximation of
%  the Bromwich integral, followed by sequence acceleration.
%  See eq. (2) in the reference below.

%                                   Andre Weideman, Updated: November 1997

%  Input:  
%      F:  transform to be inverted; string function of s,
%          e.g., F = '1./sqrt(s.^2+1)';
%      M:  number of terms emloyed in summation (50 < M < 100 is recommended.)
%      t:  value of t at which f(t) is to be computed; so far only scalar
%          t is allowed.
%     at:  free parameter (7 < at < 12 is recommended.)

% Output:  approximate value of f(t).

% Reference:   D.G. Duffy, "On the Numerical Inversion of Laplace Transforms:
%              Comparison of Three New Methods on Characteristic Problems 
%              from Applications", ACM TOMS, Vol. 19, pp. 333-359 (1993).

%  Preliminaries:

a  = at/t;
n  = 1:M;
sg = (-1).^n;                     

%  Evaluate the transform

s  = a;                      F0 = real(eval(F));
s  = a + n*pi*i/t;           Fr = real(eval(F));
s  = a + (2*n-1)*pi*i/(2*t); Fi = imag(eval(F));

%  Compute partial sums

aa = cumsum((sg.*(Fr+Fi)));

%  Apply sequence accelaration (epsilon scheme) to last 50% of partial sums.

aa = aa(length(aa)-M/2:length(aa));                
sa = wynneps(aa);                                  

%  Compute final answer:

f  = exp(at)/(2*t)*(F0/2 + sa);














