function [t,LE]=FO_NC_Lyapunov(ne,ext_fcn,t_start,t_end,h_norm, x_start,h,q,out);
% 
% Program to compute the spectrum of Lyapunov exponents (LEs) as function of time of systems of
% noncommensurate fractional order (FO) defined with Caputo's derivative 
% June 2021
% author Marius-F. Danca
% email: danca@rist.ro

% The program uses a fast variant of the Adams-Bashforth-Moulton for
% noncommensurate fractional-order differential equations: fde_pi12_pc.m , by Roberto Garrappa:
% shorturl.at/izJW3, the routine being presented in 
% R.Garrappa,  Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial, Mathematics 2018, 6(2), 16 
% m-files required: fde_pi12_pc, FO_NC_Lyapunov.m and the function containing the extended system
% Files FO_NC_Lyapunov.m, fde_pi12_pc.m and LE_FO_RF.m must be in the same folder.
% FO_NC_Lyapunov.m was developed, by modifying the program for the commensurate order 
% FO_Lyapunov.m: shorturl.at/suzF7 (see also Lyapunov.m by V.N. Govorukhin)
% 
% Input:
% ne - system dimension;
% ext_fcn - function containing the extended system of FO (see e.g. LE_RF.m in [1] or [2]);
% t_start, t_end - time span (fde_pi12_pc);
% h_norm - step for Gram-Schmidt renormalization;
% x_start - initial condition;
% h - integration step for fde_pi12_pc;
% q=[q_1;q_2;...;q_ne] - the fractional order;
% out - priniting step of LEs values;
% out==0 - no print.
%
% Output:
% t - time values;
% LE Lyapunov exponents to each time value printed every 'out' step.
%
% How to use it:
%
% [t,LE]=FO_NC_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out);
%
% Example: 
%
% [t,LE]=FO_NC_Lyapunov(3,@LE_FO_RF,0,.5,1000,[.1,.1,.1]',0.01,[.92,.92,.92],500)
%
%==================================================
% Cite the code as:
%
% Marius-F. Danca, Matlab code for Lyapunov 
% exponents of fractional-order systems, 
% Part II: The non-commensurate case, IJBC, 
% accepted on Jun 12, 2021.
%
% The algorithm for COMMENSURATE order is
% explained in:
%
% [1] Marius-F. Danca, Nikolay Kuznetsov, Matlab
% code for Lyapunov exponents of fractional order
% systems, IJBC, 28(05), 1850067 (2018)
%
% =================================================
figure;
hold on; 
% Set orders q for extended system (ne system equations + ne*ne variational
% equations
 q=repmat(q',ne+1,1);
% Memory allocation 
x=zeros(ne*(ne+1),1); 
x0=x;
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);

% Initial values for extended system
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
LE=zeros(ne,1);
% Main loop
it=1;
while it<=n_it
% Solutuion of extended ODE system of FO using FDE12 routine
    [T,Y] = fde_pi12_pc(q,ext_fcn,t,t+h_norm,x,h); 
    t=t+h_norm;
    Y=transpose(Y); 
    x=Y(size(Y,1),:); %solution at t+h_norm
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
%   construct new orthonormal basis by gram-schmidt
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)^2; 
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1); 
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)^2;
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
            j=j+1;
    end
%   update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
%   normalize exponent
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x(ne*j+i)=x0(ne*i+j);
            j=j+1;
        end
        i=i+1;
    end;
        x=transpose(x);
        it=it+1;
%       print and plot results
        if (mod(it,out)==0)
            fprintf('%10.2f %10.4f %10.4f %10.4f\n',[t,LE']);
        end;
        plot(t,LE)
end
xlabel('t','fontsize',16)
ylabel('LEs','fontsize',14)
set(gca,'fontsize',14)
box on;
line([0,t],[0,0],'color','k')
axis tight
