function F = NumericalIntg( t,dF,Fo )
%__________________________________________________________________________
% in_args: t, time vector
%          dF, function to integrate
%          Fo, initial value
% out_args: dF, first derivative
%          ddF, second derivative
%
% Specifies the differentce method 
 INTG_meth = 1; % Two-point estimation (Newton's difference quocient)
                % Three-point estimation: secant line
                % Five-point method
 Ts = t(2)-t(1);
 DIM = size(F,2);
  dF = zeros(size(F,1),DIM);     % Pre-allocating for speed
 ddF = zeros(size(F,1),DIM);     % Pre-allocating for speed
 switch INTG_meth,
     case 0,
             dF = ([F(:,2:end) F(:,end)]-F )/Ts;
            ddF = ([dF(:,2:end) dF(:,end)]-dF )/Ts; 
     case 1,            
            for kk = 1:DIM, % First derivative
                if kk == 1,
                   dF(:,kk) = ( F(:,kk+1)-F(:,kk) )/Ts;                      
                elseif kk == DIM
                       dF(:,kk) = ( F(:,kk)-F(:,kk-1) )/Ts;
                else
                    dF(:,kk) = ( F(:,kk+1)-F(:,kk-1) )/(2*Ts);
                end                 
            end
            for kk = 1:DIM, % Second derivative
                if kk == 1,
                   ddF(:,kk) = ( dF(:,kk+1)-dF(:,kk) )/Ts;                      
                elseif kk == DIM
                       ddF(:,kk) = ( dF(:,kk)-dF(:,kk-1) )/Ts;
                else
                    ddF(:,kk) = ( dF(:,kk+1)-dF(:,kk-1) )/(2*Ts);
                end                 
            end
     case 2,
            for kk = 1:DIM, % First derivative
                if kk == 1,
                   dF(:,kk) = ( F(:,kk+1)-F(:,kk) )/Ts;
                elseif kk == 2 | kk == 3 | kk == DIM-1 | kk == DIM-2,
                       dF(:,kk) = ( F(:,kk+1)-F(:,kk-1) )/(2*Ts);
                elseif kk == DIM
                       dF(:,kk) = ( F(:,kk)-F(:,kk-1) )/Ts;
                else
                    dF(:,kk) = ( -F(:,kk+2)+8*F(:,kk+1)-8*F(:,kk-1)+F(:,kk-2) )/(12*Ts);
                end                 
            end
             for kk = 1:DIM, % Second derivative
                if kk == 1,
                   ddF(:,kk) = ( dF(:,kk+1)-dF(:,kk) )/Ts;
                elseif kk == 2 | kk == 3 | kk == DIM-1 | kk == DIM-2,
                       ddF(:,kk) = ( dF(:,kk+1)-dF(:,kk-1) )/(2*Ts);
                elseif kk == DIM
                       ddF(:,kk) = ( dF(:,kk)-dF(:,kk-1) )/Ts;
                else
                    ddF(:,kk) = ( -dF(:,kk+2)+8*dF(:,kk+1)-8*dF(:,kk-1)+dF(:,kk-2) )/(12*Ts);
                end                 
            end
            
 end
%__________________________________________________________________________
end

