stefhest.m

function s = stefhest(F,x,t,N)
%STEFHEST - Numerical Laplace inversion with the Stefhest method
%
% Syntax: s = stefhest(F,x,t,N)
%
%   F = name of a function in the laplace domain
%       F should be of the form F(x,p) = L( F(x,t) )
%   x = vector of the parameters of the function
%   t = time vector 
%   N = parameter of the Stefhest algorithm (must be even)
%   N = 12 by default
%
% Description:
%   Numerical Laplace inversion with the Stefhest method.
%
% Example:
%   sd = stefhest('theis_l',[],td)
%   sd = stefhest('cooper_l',[1,1e-2],td,8)
%
% See also: dehoog
%



if(nargin==3)           % Default value for N if not given by user
  N=12;
end

% Check the time vector

if size(t,1) == 1        % It is a line vector then no problem 
  flag=0;
elseif size(t,2) == 1    % It is a row vector then we have to transpose it 
  flag=1;
  t=t';
else                     % It is a matrix, it cannot be treated
  error('Time Vector cannot be a matrix')
end


for i=1:N               % Computation of the coefficient Vi
   
   vi_dummy = 0;
   for k=fix((i+1)/2):min(i,N/2) 
      
      vi_dummy = vi_dummy + (k.^(N/2)*prod(1:2*k))./...
         (prod(1:N/2-k).*prod(1:k).*prod(1:k-1).*...
         prod(1:i-k).*prod(1:2*k-i));
   end
   
   V(i) = (-1).^((N/2)+i)*vi_dummy;
   
end

i=1:N;                  % Vecteur ligne des valeurs de 1 a N

p=(i*log(2))'*(1./t);   % Matrice 
                        % Chaque colonne correspond a un temps
                        %  (provient du vecteur ligne t)
                        % Chaque ligne a une valeur de i
                        % p(i,j) = i*log(2)/t(j)



up=feval(F,x,p);        % Calcul de la function dans le domaine 
                        % de Laplace pour chaque p(i,j)

                        % up est une matrice de meme structure 
                        % que p 
                        % up(i,j) = F( p(i,j) )


%for k=1:N               % Multiplication de chaque ligne de up
%  su(k,:)=V(k).*up(k,:);% par le coefficient Vi
%end

VV=repmat(V',1,size(up,2)); % Alternative equivalent of the previous
su=VV.*up;                  % equation in matrix format

s=log(2)./t.*sum(su);

if flag==1               % We have to transpose the result so that
  s=s';                  % it can be plotted together with t
end