function f = invlap(F, tt, alpha);
% f = invlap(F, t, alpha);
% alpha   largest pole of F (default zero)
% f       vector of real-space values f(t)
if nargin <= 2, alpha = 0; end
tol = 1e-9; f = []; t = tt(:);
allt = t;	
logallt = log10(allt);
iminlogallt = floor(min(logallt));
imaxlogallt = ceil(max(logallt));
for ilogt = iminlogallt:imaxlogallt,	 
    t = allt(find((logallt>=ilogt) & (logallt<(ilogt+1))));
    if ~isempty(t),	
        T = max(t)*2;
%----------------------------------------
% modified by B. Welfert
        gamma = alpha-log(tol)/(2*T);
%        gamma = alpha-1.2*log(tol)/(2*T);
%----------------------------------------
        nt = length(t);
        M = 20;
        run = [0:1:2*M]'; 
        s = gamma + i*pi*run/T;
        a = feval(F,s);
        a(1) = a(1)/2; 
        e = zeros(2*M+1, M+1);
        q = zeros(2*M  , M+1); 
        e(:,1) = zeros(2*M+1,1);
        q(:,2) = a(2:2*M+1,1)./a(1:2*M,1);
        for r = 2:M+1, 
            e(1:2*(M-r+1)+1,r) = ...
                q(2:2*(M-r+1)+2,r) - q(1:2*(M-r+1)+1,r) + e(2:2*(M-r+1)+2,r-1);
            if r<M+1,    
                rq = r+1;
                q(1:2*(M-rq+1)+2,rq) = ...
                    q(2:2*(M-rq+1)+3,rq-1).*e(2:2*(M-rq+1)+3,rq-1)./e(1:2*(M-rq+1)+2,rq-1);
            end
        end
        d = zeros(2*M+1,1);
        d(1,1) = a(1,1);
        d(2:2:2*M,1) = -q(1,2:M+1).'; 
        d(3:2:2*M+1,1) = -e(1,2:M+1).'; 
        A = zeros(2*M+2,nt);
        B = zeros(2*M+2,nt);
        A(2,:) = d(1,1)*ones(1,nt);
        B(1:2,:) = ones(2,nt);
        z = exp(i*pi*t'/T);
        for n = 3:2*M+2,
            A(n,:) = A(n-1,:) + d(n-1,1)*ones(1,nt).*z.*A(n-2,:); 
            B(n,:) = B(n-1,:) + d(n-1,1)*ones(1,nt).*z.*B(n-2,:); 
        end
        h2M = .5 * ( ones(1,nt) + ( d(2*M,1)-d(2*M+1,1) )*ones(1,nt).*z );
        R2Mz = -h2M.*(ones(1,nt) - ...
            (ones(1,nt)+d(2*M+1,1)*ones(1,nt).*z/(h2M).^2).^.5);
        A(2*M+2,:) = A(2*M+1,:) + R2Mz .* A(2*M,:);
        B(2*M+2,:) = B(2*M+1,:) + R2Mz .* B(2*M,:);
        fpiece = ( 1/T * exp(gamma*t') .* real(A(2*M+2,:)./B(2*M+2,:)) )';
        f = [f; fpiece];
    end
end
if size(tt,1)==1, f=f'; end