Вычисление Ляпуновских показателей

clear;
clc;
 
 
a = 1.4;
b = -0.3;
summa = zeros(2,1);
n = 600;
 
 
x(1) = 0.001; %начальные условия для самого отображения
y(1) = 0.001;
 
xx(1) = 1; %начальные условия для первой системы уравнений в вариациях
yy(1) = 0;
 
X(1) = 0; %начальные условия для второй системы уравнений в вариациях
Y(1) = 1;
 
    for j=1:n
        x(j+1) = 1-a*x(j)^2-b*y(j); %решаем численно само отображение 
        y(j+1) = x(j);
        
        xx(j+1) = -2*a*x(j)*xx(j)-b*yy(j); %решаем численно первую систему
        yy(j+1) = xx(j);
        
        X(j+1) = -2*a*x(j)*X(j)-b*Y(j); %решаем численно вторую систему
        Y(j+1) = X(j);
        
        k = [ xx(j+1), yy(j+1) ]; %эти два вектора создаю для удобства
        l = [ X(j+1), Y(j+1) ];
        
        w = k; %ортогонализация векторов
        z = l - dot(l, w)*w;
        
        wNorm = norm( w ); %соостветсвенно нормировка векторов
        zNorm = norm( z );
        
        summa(1) = summa(1) + log( wNorm ); %вычисляем сумму для подсчета показателей Ляпунова
        summa(2) = summa(2) + log( zNorm );
        
        xx(j+1) = w(1)/wNorm; %нормируем вектора для следующей итерации
        yy(j+1) = w(2)/wNorm;
        
        X(j+1) = z(1)/zNorm; %нормируем вектора для следующей итерации
        Y(j+1) = z(2)/zNorm;
    end;
 
lamda1 = summa(1)/n; %собственно сами показателя Ляпунова
lamda2 = summa(2)/n;