Problems computing the Lyapunov spectrum depending on time interval T

Question

I intent to determine the Lyapunov spectrum of nonlinear dynamic systems. To do so, I implemented some code in MATLAB, based on algorithms from literature and very similar to the way, described in this thread: Numerical calculation of Lyapunov exponents without Jacobian

I added my own MATLAB code at the end of this article. My problem now is, that the results depends very much on the length you choose for the time interval T, or better said tspan, you realize the integration for, in every iteration. In the literature a proper length of the interval is estimated as 10 to twenty times the natural period of the system. But this does not work for me. In most cases this interval already seems to be to big and I obtain wrong Lyapunov exponents. In some cases it works with smaller intervals but I can not see a pattern there, sometimes it seems to be impossible to find the right interval and the Lyapunov exponents change with every variation of the time interval.

Therefore I hoped to find someone here, you might have experienced a similar problem or has an idea where this strong dependence of the time interval comes from and how to handle it, to obtain the right Lyapunov Exponents.

Thank you very much for your help already!

function lyapunov = lyapunov_spektrum(dgl, tspan, x0, perturbation, kmax, tolr, tola)

% set variables and parameters
k=0;                                % set counter for iterations
dim = length(x0);                   % get dimension from length of x0
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
I = eye(dim, dim);
dy = perturbation;                  
y = zeros(dim, dim);                % matrix with vectors of perturbations

for zaehler_0 = 1:dim
    x0_d(:,zaehler_0) = x0;         % building perturbated initial conditions
end;
x0_d = x0_d - I*dy;

lambda = zeros(length(x0),1);       % saves current lyapunov exponent in the end of every iteration
sum = zeros(length(x0),1);          


% begin of iterations to compute lyapunov spectrum

while k < kmax
    k = k+1;
    lambda_alt = lambda;
    
    % computing solution to initial coniditon x0

    [t, xsol] = ode15s(str2func(dgl), tspan, x0,opts);
    x = xsol(end,:);                            % Vektor mit Lösung am Intervallende
    x0 = x;                                     % Übergabe der AW für nächste k-te Iteration
    
    
    % computing solutions to every single perutbated initial condition (cloumns of x0_d)
    
    for zaehler_a = 1:dim
        [t, xsol_d] = ode15s(str2func(dgl), tspan, x0_d(:,zaehler_a), opts);   
        x_d(:,zaehler_a) = xsol_d(end,:)';      % save value at the end of time interval tspan
    end;
    
    
    % building perturbation vectors (columns of y)
    
    for zaehler_b = 1:dim
        y(:, zaehler_b) = x' - x_d(:,zaehler_b);                               
    end;
    
    
   
    % Gram-Schmidt procedure to orthonormalize perturbation vectors
    
    for zaehler_3 = 1:dim
        v(:,zaehler_3) = y(:,zaehler_3);                                       
        for zaehler_4 = 1:zaehler_3 - 1
            v(:,zaehler_3) = v(:,zaehler_3) - dot(v(:,zaehler_3), u(:,zaehler_4))*u(:,zaehler_4);
        end;
        u(:,zaehler_3) = (v(:,zaehler_3)/norm(v(:,zaehler_3)));             
        sum(zaehler_3) = sum(zaehler_3) + log2(1/dy*norm(v(:,zaehler_3)));
        lambda(zaehler_3) = 1/(k*tspan(end))*sum(zaehler_3);   % calculating lyapunov exponents
    end;
    
   
    % building new initial conditions for the next iteration

    for zaehler_5 = 1:dim                                      
        x0_d(:,zaehler_5) = x0' - dy*u(:,zaehler_5);                              
        
        
    % termination criterion
    
    if norm(lambda_alt-lambda) < tolr*norm(lambda) + tola       
        lyapunov = lambda;                                                      
        disp('The lyapunov exponents are');
        return;
    end;
end;


% if there are no lyapunov exponents found

disp('There were no lyapunov exponents found.');