%5.3  ilaplacecode.m

%This code was contributed by Andre Weideman

function f = talbotmod(F,t,N)

%   The function f = talbot(F,t,N) implements Talbot's
%   *modified* method for the inversion of the Laplace 
%   Transform.   F is the transform to be inverted, supplied
%   as an inline function of z.   (The code assumes that F is 
%   real when z is real (*).)   N is the number of F evaluations.   
%   The method is based on a 2N panel midpoint rule approximation, 
%   but because of the symmetry (*) it can be reduced to N evaluations 
%   of F.   The parameters sigma, mu, and nu are tuned for optimality
%   in the case where all the singularities of F are on or near the 
%   negative imaginary axis.   In this case a convergence rate of 
%   about exp(-2.7N) can be expected, i.e.,  N = 14 should usually be 
%   sufficient for full precision in MATLAB.  The code does not yet 
%   make provision for multiple values of t. 
%
%   Example of usage:
%    >> F = inline('1./(z.*sqrt(z+1))')
%    >> f = talbotmod(F,1,10) 
%                                      JAC Weideman, September 2006

sigma = -1.2244*N/t; mu = 1.0035*N/t; nu = 0.5272;  % define optimal params
alpha = 0.6407;                                     % new parameter
theta = (2*[0:N-1]+1)*pi/(2*N);                     % equispaced theta vals
z  = sigma + mu*(theta.*cot(alpha*theta) + nu*i*theta);   % modified contour
dz = mu*(cot(alpha*theta)-alpha*theta.*csc(alpha*theta).^2 + nu*i); % dz/dtheta
f  =  1/N*imag(sum(exp(z*t).*F(z).*dz));            % midpoint sum
 
y