[ft,t]=nilt(F,tm);

% NILT numerical inverse Laplace transform
% except from Programs for Fast Numerical Inversion of Laplace 
% Transforms in Matlab Language Environment, Lubomir Brancik, 
% Konference MATLAB 99 ZCU, Plzen, 1999, pp. 27-39
%
% function [ft, t] = nilt(F, tm)
% 
% F: file name of the transfer function
% tm: time range in which to calculate the inverse Laplace transform
% ft: vector of the nemerical value of the inverse Laplace transform
% t: vector of the time, t(1) = 0, t(end) = tm

function [ft,t]=nilt(F,tm);

%alfa=0; M=256; P=2;
%alfa=0; M=512; P=2;
alfa=0; M=1024; P=2;
%alfa=0; M=2048; P=2;

% In boundary control problems, if the solution includes very high
% frequency components, such as boundary control of beam equation with
% time delays in the tip velocity measurement using static controller 
% only, experience shows M should be at least 2048. Otherwise M = 1024
% should be enough. The values of alfa and P are not tuned so far. 
% Commented by Jinsong Liang


N=2*M; wyn=2*P+1;
t=linspace(0,tm,M);

NT=2*tm*N/(N-2); omega=2*pi/NT;
c=alfa+25/NT; s=c-i*omega*(0:N+wyn-2);
Fsc=feval(F,s);
ft=fft(Fsc(1:N)); ft=ft(1:M);
for n=N:N+wyn-2
ft(n-N+2,:)=Fsc(n+1)*exp(-i*n*omega*t);
end
ft1=cumsum(ft); ft2=zeros(wyn-1,M);
for I=1:wyn-2
ft=ft2+1./diff(ft1);
ft2=ft1(2:wyn-I,:); ft1=ft;
end
ft=ft2+1./diff(ft1); ft=2*real(ft)-Fsc(1);
ft=exp(c*t)/NT.*ft; ft(1)=2*ft(1);
%plot(t,ft); xlabel('t'); ylabel('f(t)'); grid on;