# Lyapunov exponents matlabself.__wrap_n=self.__wrap_n||(self.CSS&&CSS.supports("text-wrap","balance")?1:2);self.__wrap_b=(e,t,r)=>{let o=(r=r||document.querySelector(`[data-br="\${e}"]`)).parentElement,n=e=>r.style.maxWidth=e+"px";r.style.maxWidth="";let a=o.clientWidth,i=o.clientHeight,s=a/2-.25,l=a+.5,u;if(a){for(n(s),s=Math.max(r.scrollWidth,s);s+1<l;)n(u=Math.round((s+l)/2)),o.clientHeight===i?l=u:s=u;n(l*t+a*(1-t))}r.__wrap_o||"undefined"!=typeof ResizeObserver&&(r.__wrap_o=new ResizeObserver(()=>{self.__wrap_b(0,+r.dataset.brr,r)})).observe(o)};self.__wrap_n!=1&&self.__wrap_b(":R1emmlkqla:",1)

... views fde12.m
``````function [t, y] = fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)
%FDE12  Solves an initial value problem for a non-linear differential
%       equation of fractional order (FDE). The code implements the
%       predictor-corrector PECE method of Adams-Bashforth-Moulton type
%       described in .
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value
%   problem for the FDE, or the system of FDEs, of order ALPHA > 0
%      D^ALPHA Y(t) = FDEFUN(T,Y(T))
%      Y^(k)(T0) = Y0(:,k+1), k=0,...,m-1
%   where m is the smallest integer grater than ALPHA and D^ALPHA is the
%   fractional derivative according to the Caputo's definition. FDEFUN is a
%   function handle corresponding to the vector field of the FDE and for a
%   scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The
%   set of initial conditions Y0 is a matrix with a number of rows equal to
%   the size of the problem (hence equal to the number of rows of the
%   output of FDEFUN) and a number of columns depending on ALPHA and given
%   by m. The step-size H>0 is assumed constant throughout the integration.
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with
%   the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with
%   the selected number MU of multiple corrector iterations. The following
%   values for MU are admissible:
%   MU = 0 : the corrector is not evaluated and the solution is provided
%   just by the predictor method (the first order rectangular rule);
%   MU > 0 : the corrector is evaluated by the selected number MU of times;
%   the classical PECE method is obtained for MU=1;
%   MU = Inf : the corrector is evaluated for a certain number of times
%   until convergence of the iterations is reached (for convergence the
%   difference between two consecutive iterates is tested).
%   The defalut value for MU is 1
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to
%   specify the tolerance for testing convergence when MU = Inf. If not
%   specified, the default value MU_TOL = 1.E-6 is used.
%
%   FDE12 is an implementation of the predictor-corrector method of
%   Adams-Bashforth-Moulton studied in . Convergence and accuracy of
%   the method are studied in . The implementation with multiple
%   corrector iterations has been proposed and discussed for multiterm FDEs
%   in . In this implementation the discrete convolutions are evaluated
%   by means of the FFT algorithm described in  allowing to keep the
%   computational cost proportional to N*log(N)^2 instead of N^2 as in the
%   classical implementation; N is the number of time-point in which the
%   solution is evaluated, i.e. N = (TFINAL-T)/H. The stability properties
%   of the method implemented by FDE12 have been studied in .
%
%    K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical
%   solution of differential equations of fractional order, in: S. Heinzel,
%   T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998,
%   Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999,
%   pp. 57-71.
%
%    K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a
%   fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.
%
%    K. Diethelm, Efficient solution of multi-term fractional
%   differential equations using P(EC)mE methods, Computing 71 (2003), pp.
%   305-319.
%
%    E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of
%   nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput.
%   6 (3) (1985) 532-541.
%
%    R. Garrappa, On linear stability of predictor-corrector algorithms
%   for fractional differential equations, Internat. J. Comput. Math. 87
%   (10) (2010) 2281-2290.
%
%   Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy
%   garrappa at dm dot uniba dot it
%   Revision: 1.2 - Date: July, 6 2012

% Check inputs
if nargin < 9
mu_tol = 1.0e-6 ;
if nargin < 8
mu = 1 ;
if nargin < 7
param = [] ;
end
end
end

% Check order of the FDE
if alpha <= 0
error('MATLAB:fde12:NegativeOrder',...
['The order ALPHA of the FDE must be positive. The value ' ...
'ALPHA = %f can not be accepted. See FDE12.'], alpha);
end

% Check the step--size of the method
if h <= 0
error('MATLAB:fde12:NegativeStepSize',...
['The step-size H for the method must be positive. The value ' ...
'H = %e can not be accepted. See FDE12.'], h);
end

% Structure for storing initial conditions
ic.t0 = t0 ;
ic.y0 = y0 ;
ic.m_alpha = ceil(alpha) ;
ic.m_alpha_factorial = factorial(0:ic.m_alpha-1) ;

% Structure for storing information on the problem
Probl.ic = ic ;
Probl.fdefun = fdefun ;
Probl.problem_size = size(y0,1) ;
Probl.param = param ;

% Check number of initial conditions
if size(y0,2) < ic.m_alpha
error('MATLAB:fde12:NotEnoughInputs', ...
['A not sufficient number of assigned initial conditions. ' ...
'Order ALPHA = %f requires %d initial conditions. See FDE12.'], ...
alpha,ic.m_alpha);
end

% Check compatibility size of the problem with size of the vector field
f_temp = f_vectorfield(t0,y0(:,1),Probl) ;
if Probl.problem_size ~= size(f_temp,1)
error('MATLAB:fde12:SizeNotCompatible', ...
['Size %d of the problem as obtained from initial conditions ' ...
'(i.e. the number of rows of Y0) not compatible with the ' ...
'size %d of the output of the vector field FDEFUN. ' ...
'See FDE12.'], Probl.problem_size,size(f_temp,1));
end

% Number of points in which to evaluate weights and solution
r = 16 ;
N = ceil((tfinal-t0)/h) ;
Nr = ceil((N+1)/r)*r ;
Q = ceil(log2(Nr/r)) - 1 ;
NNr = 2^(Q+1)*r ;

% Preallocation of some variables
y = zeros(Probl.problem_size,N+1) ;
fy = zeros(Probl.problem_size,N+1) ;
zn_pred = zeros(Probl.problem_size,NNr+1) ;
if mu > 0
zn_corr = zeros(Probl.problem_size,NNr+1) ;
else
zn_corr = 0 ;
end

% Evaluation of coefficients of the PECE method
nvett = 0 : NNr+1 ;
nalpha = nvett.^alpha ;
nalpha1 = nalpha.*nvett ;
PC.bn = nalpha(2:end) - nalpha(1:end-1) ;
PC.an = [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) + nalpha1(3:end)) ] ;
PC.a0 = [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;
PC.halpha1 = h^alpha/gamma(alpha+1) ;
PC.halpha2 = h^alpha/gamma(alpha+2) ;
PC.mu = mu ; PC.mu_tol = mu_tol ;

% Initializing solution and proces of computation
t = t0 + (0 : N)*h ;
y(:,1) = y0(:,1) ;
fy(:,1) = f_temp ;
[y, fy] = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;

% Main process of computation by means of the FFT algorithm
ff = [0 2 ] ; nx0 = 0 ; ny0 = 0 ;
for q = 0 : Q
L = 2^q ;
[y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, ...
zn_pred, zn_corr, N, PC, Probl ) ;
ff = [ff ff]  ; ff(end) = 4*L ;
end

% Evaluation solution in TFINAL when TFINAL is not in the mesh
if tfinal < t(N+1)
c = (tfinal - t(N))/h ;
t(N+1) = tfinal ;
y(:,N+1) = (1-c)*y(:,N) + c*y(:,N+1) ;
end
t = t(1:N+1) ; y = y(:,1:N+1) ;

end

% =========================================================================
% =========================================================================
% r : dimension of the basic square or triangle
% L : factor of resizing of the squares
function [y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...
zn_pred, zn_corr, N , PC, Probl)

nxi = nx0 ; nxf = nx0 + L*r - 1 ;
nyi = ny0 ; nyf = ny0 + L*r - 1 ;
is = 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;

i_triangolo = 0 ; stop = 0 ;
while ~stop

stop = nxi+r-1 == nx0+L*r-1 | (nxi+r-1>=Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato

[zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;

[y, fy] = Triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;
i_triangolo = i_triangolo + 1 ;

if ~stop
if nxi+r-1 == nxf   % Il triangolo finisce dove finisce il quadrato -> si scende di livello
i_Delta = ff(i_triangolo) ;
Delta = i_Delta*r ;
nxi = s_nxf(is)+1 ; nxf = s_nxf(is)  + Delta ;
nyi = s_nxf(is) - Delta +1; nyf = s_nxf(is)  ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
else % Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r  ;
is = is + 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
end
end

end

end

% =========================================================================
% =========================================================================
function [zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)

coef_beg = nxi-nyf ; coef_end = nxf-nyi+1 ;
funz_beg = nyi+1 ; funz_end = nyf+1 ;

% Evaluation convolution segment for the predictor
vett_coef = PC.bn(coef_beg:coef_end) ;
vett_funz = [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
zzn_pred = real(FastConv(vett_coef,vett_funz)) ;
zn_pred(:,nxi+1:nxf+1) = zn_pred(:,nxi+1:nxf+1) + zzn_pred(:,nxf-nyf+1-1:end-1) ;

% Evaluation convolution segment for the corrector
if PC.mu > 0
vett_coef = PC.an(coef_beg:coef_end) ;
if nyi == 0 % Evaluation of the lowest square
vett_funz = [zeros(Probl.problem_size,1) , fy(:,funz_beg+1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
else % Evaluation of any square but not the lowest
vett_funz = [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
end
zzn_corr = real(FastConv(vett_coef,vett_funz)) ;
zn_corr(:,nxi+1:nxf+1) = zn_corr(:,nxi+1:nxf+1) + zzn_corr(:,nxf-nyf+1:end) ;
else
zn_corr = 0 ;
end

end

% =========================================================================
% =========================================================================
function [y, fy] = Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)

for n = nxi : min(N,nxf)

% Evaluation of the predictor
Phi = zeros(Probl.problem_size,1) ;
if nxi == 1 % Case of the first triangle
j_beg = 0 ;
else % Case of any triangle but not the first
j_beg = nxi ;
end
for j = j_beg : n-1
Phi = Phi + PC.bn(n-j)*fy(:,j+1) ;
end

St = StartingTerm(t(n+1),Probl.ic) ;
y_pred = St + PC.halpha1*(zn_pred(:,n+1) + Phi) ;
f_pred = f_vectorfield(t(n+1),y_pred,Probl) ;

if PC.mu == 0
y(:,n+1) = y_pred ;
fy(:,n+1) = f_pred ;
else
j_beg = nxi ;
Phi = zeros(Probl.problem_size,1) ;
for j = j_beg : n-1
Phi = Phi + PC.an(n-j+1)*fy(:,j+1) ;
end
Phi_n = St + PC.halpha2*(PC.a0(n+1)*fy(:,1) + zn_corr(:,n+1) + Phi) ;
yn0 = y_pred ; fn0 = f_pred ;
stop = 0 ; mu_it = 0 ;
while ~stop
yn1 = Phi_n + PC.halpha2*fn0 ;
mu_it = mu_it + 1 ;
if PC.mu == Inf
stop = norm(yn1-yn0,inf) < PC.mu_tol ;
if mu_it > 100 && ~stop
warning('MATLAB:fde12:NonConvegence',...
strcat('It has been requested to run corrector iterations until convergence but ', ...
'the process does not converge to the tolerance %e in 100 iteration'),PC.mu_tol) ;
stop = 1 ;
end
else
stop = mu_it == PC.mu ;
end
fn1 = f_vectorfield(t(n+1),yn1,Probl) ;
yn0 = yn1 ; fn0 = fn1 ;
end
y(:,n+1) = yn1 ;
fy(:,n+1) = fn1 ;
end
end

end

% =========================================================================
% =========================================================================
function z = FastConv(x, y)

r = length(x) ; problem_size = size(y,1) ;

z = zeros(problem_size,r) ;
X = fft(x,r) ;
for i = 1 : problem_size
Y = fft(y(i,:),r) ;
Z = X.*Y ;
z(i,:) = ifft(Z,r) ;
end

end

% =========================================================================
% =========================================================================
function f = f_vectorfield(t,y,Probl)

if isempty(Probl.param)
f = feval(Probl.fdefun,t,y) ;
else
f = feval(Probl.fdefun,t,y,Probl.param) ;
end

end

% =========================================================================
% =========================================================================
function ys = StartingTerm(t,ic)

ys = zeros(size(ic.y0,1),1) ;
for k = 1 : ic.m_alpha
ys = ys + (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;
end

end``````
fde12.m
``````function [t, y] = fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)
%FDE12  Solves an initial value problem for a non-linear differential
%       equation of fractional order (FDE). The code implements the
%       predictor-corrector PECE method of Adams-Bashforth-Moulton type
%       described in .
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value
%   problem for the FDE, or the system of FDEs, of order ALPHA > 0
%      D^ALPHA Y(t) = FDEFUN(T,Y(T))
%      Y^(k)(T0) = Y0(:,k+1), k=0,...,m-1
%   where m is the smallest integer grater than ALPHA and D^ALPHA is the
%   fractional derivative according to the Caputo's definition. FDEFUN is a
%   function handle corresponding to the vector field of the FDE and for a
%   scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The
%   set of initial conditions Y0 is a matrix with a number of rows equal to
%   the size of the problem (hence equal to the number of rows of the
%   output of FDEFUN) and a number of columns depending on ALPHA and given
%   by m. The step-size H>0 is assumed constant throughout the integration.
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with
%   the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with
%   the selected number MU of multiple corrector iterations. The following
%   values for MU are admissible:
%   MU = 0 : the corrector is not evaluated and the solution is provided
%   just by the predictor method (the first order rectangular rule);
%   MU > 0 : the corrector is evaluated by the selected number MU of times;
%   the classical PECE method is obtained for MU=1;
%   MU = Inf : the corrector is evaluated for a certain number of times
%   until convergence of the iterations is reached (for convergence the
%   difference between two consecutive iterates is tested).
%   The defalut value for MU is 1
%
%   [T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to
%   specify the tolerance for testing convergence when MU = Inf. If not
%   specified, the default value MU_TOL = 1.E-6 is used.
%
%   FDE12 is an implementation of the predictor-corrector method of
%   Adams-Bashforth-Moulton studied in . Convergence and accuracy of
%   the method are studied in . The implementation with multiple
%   corrector iterations has been proposed and discussed for multiterm FDEs
%   in . In this implementation the discrete convolutions are evaluated
%   by means of the FFT algorithm described in  allowing to keep the
%   computational cost proportional to N*log(N)^2 instead of N^2 as in the
%   classical implementation; N is the number of time-point in which the
%   solution is evaluated, i.e. N = (TFINAL-T)/H. The stability properties
%   of the method implemented by FDE12 have been studied in .
%
%    K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical
%   solution of differential equations of fractional order, in: S. Heinzel,
%   T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998,
%   Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999,
%   pp. 57-71.
%
%    K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a
%   fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.
%
%    K. Diethelm, Efficient solution of multi-term fractional
%   differential equations using P(EC)mE methods, Computing 71 (2003), pp.
%   305-319.
%
%    E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of
%   nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput.
%   6 (3) (1985) 532-541.
%
%    R. Garrappa, On linear stability of predictor-corrector algorithms
%   for fractional differential equations, Internat. J. Comput. Math. 87
%   (10) (2010) 2281-2290.
%
%   Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy
%   garrappa at dm dot uniba dot it
%   Revision: 1.2 - Date: July, 6 2012

% Check inputs
if nargin < 9
mu_tol = 1.0e-6 ;
if nargin < 8
mu = 1 ;
if nargin < 7
param = [] ;
end
end
end

% Check order of the FDE
if alpha <= 0
error('MATLAB:fde12:NegativeOrder',...
['The order ALPHA of the FDE must be positive. The value ' ...
'ALPHA = %f can not be accepted. See FDE12.'], alpha);
end

% Check the step--size of the method
if h <= 0
error('MATLAB:fde12:NegativeStepSize',...
['The step-size H for the method must be positive. The value ' ...
'H = %e can not be accepted. See FDE12.'], h);
end

% Structure for storing initial conditions
ic.t0 = t0 ;
ic.y0 = y0 ;
ic.m_alpha = ceil(alpha) ;
ic.m_alpha_factorial = factorial(0:ic.m_alpha-1) ;

% Structure for storing information on the problem
Probl.ic = ic ;
Probl.fdefun = fdefun ;
Probl.problem_size = size(y0,1) ;
Probl.param = param ;

% Check number of initial conditions
if size(y0,2) < ic.m_alpha
error('MATLAB:fde12:NotEnoughInputs', ...
['A not sufficient number of assigned initial conditions. ' ...
'Order ALPHA = %f requires %d initial conditions. See FDE12.'], ...
alpha,ic.m_alpha);
end

% Check compatibility size of the problem with size of the vector field
f_temp = f_vectorfield(t0,y0(:,1),Probl) ;
if Probl.problem_size ~= size(f_temp,1)
error('MATLAB:fde12:SizeNotCompatible', ...
['Size %d of the problem as obtained from initial conditions ' ...
'(i.e. the number of rows of Y0) not compatible with the ' ...
'size %d of the output of the vector field FDEFUN. ' ...
'See FDE12.'], Probl.problem_size,size(f_temp,1));
end

% Number of points in which to evaluate weights and solution
r = 16 ;
N = ceil((tfinal-t0)/h) ;
Nr = ceil((N+1)/r)*r ;
Q = ceil(log2(Nr/r)) - 1 ;
NNr = 2^(Q+1)*r ;

% Preallocation of some variables
y = zeros(Probl.problem_size,N+1) ;
fy = zeros(Probl.problem_size,N+1) ;
zn_pred = zeros(Probl.problem_size,NNr+1) ;
if mu > 0
zn_corr = zeros(Probl.problem_size,NNr+1) ;
else
zn_corr = 0 ;
end

% Evaluation of coefficients of the PECE method
nvett = 0 : NNr+1 ;
nalpha = nvett.^alpha ;
nalpha1 = nalpha.*nvett ;
PC.bn = nalpha(2:end) - nalpha(1:end-1) ;
PC.an = [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) + nalpha1(3:end)) ] ;
PC.a0 = [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;
PC.halpha1 = h^alpha/gamma(alpha+1) ;
PC.halpha2 = h^alpha/gamma(alpha+2) ;
PC.mu = mu ; PC.mu_tol = mu_tol ;

% Initializing solution and proces of computation
t = t0 + (0 : N)*h ;
y(:,1) = y0(:,1) ;
fy(:,1) = f_temp ;
[y, fy] = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;

% Main process of computation by means of the FFT algorithm
ff = [0 2 ] ; nx0 = 0 ; ny0 = 0 ;
for q = 0 : Q
L = 2^q ;
[y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, ...
zn_pred, zn_corr, N, PC, Probl ) ;
ff = [ff ff]  ; ff(end) = 4*L ;
end

% Evaluation solution in TFINAL when TFINAL is not in the mesh
if tfinal < t(N+1)
c = (tfinal - t(N))/h ;
t(N+1) = tfinal ;
y(:,N+1) = (1-c)*y(:,N) + c*y(:,N+1) ;
end
t = t(1:N+1) ; y = y(:,1:N+1) ;

end

% =========================================================================
% =========================================================================
% r : dimension of the basic square or triangle
% L : factor of resizing of the squares
function [y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...
zn_pred, zn_corr, N , PC, Probl)

nxi = nx0 ; nxf = nx0 + L*r - 1 ;
nyi = ny0 ; nyf = ny0 + L*r - 1 ;
is = 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;

i_triangolo = 0 ; stop = 0 ;
while ~stop

stop = nxi+r-1 == nx0+L*r-1 | (nxi+r-1>=Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato

[zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;

[y, fy] = Triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;
i_triangolo = i_triangolo + 1 ;

if ~stop
if nxi+r-1 == nxf   % Il triangolo finisce dove finisce il quadrato -> si scende di livello
i_Delta = ff(i_triangolo) ;
Delta = i_Delta*r ;
nxi = s_nxf(is)+1 ; nxf = s_nxf(is)  + Delta ;
nyi = s_nxf(is) - Delta +1; nyf = s_nxf(is)  ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
else % Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r  ;
is = is + 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
end
end

end

end

% =========================================================================
% =========================================================================
function [zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)

coef_beg = nxi-nyf ; coef_end = nxf-nyi+1 ;
funz_beg = nyi+1 ; funz_end = nyf+1 ;

% Evaluation convolution segment for the predictor
vett_coef = PC.bn(coef_beg:coef_end) ;
vett_funz = [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
zzn_pred = real(FastConv(vett_coef,vett_funz)) ;
zn_pred(:,nxi+1:nxf+1) = zn_pred(:,nxi+1:nxf+1) + zzn_pred(:,nxf-nyf+1-1:end-1) ;

% Evaluation convolution segment for the corrector
if PC.mu > 0
vett_coef = PC.an(coef_beg:coef_end) ;
if nyi == 0 % Evaluation of the lowest square
vett_funz = [zeros(Probl.problem_size,1) , fy(:,funz_beg+1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
else % Evaluation of any square but not the lowest
vett_funz = [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
end
zzn_corr = real(FastConv(vett_coef,vett_funz)) ;
zn_corr(:,nxi+1:nxf+1) = zn_corr(:,nxi+1:nxf+1) + zzn_corr(:,nxf-nyf+1:end) ;
else
zn_corr = 0 ;
end

end

% =========================================================================
% =========================================================================
function [y, fy] = Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)

for n = nxi : min(N,nxf)

% Evaluation of the predictor
Phi = zeros(Probl.problem_size,1) ;
if nxi == 1 % Case of the first triangle
j_beg = 0 ;
else % Case of any triangle but not the first
j_beg = nxi ;
end
for j = j_beg : n-1
Phi = Phi + PC.bn(n-j)*fy(:,j+1) ;
end

St = StartingTerm(t(n+1),Probl.ic) ;
y_pred = St + PC.halpha1*(zn_pred(:,n+1) + Phi) ;
f_pred = f_vectorfield(t(n+1),y_pred,Probl) ;

if PC.mu == 0
y(:,n+1) = y_pred ;
fy(:,n+1) = f_pred ;
else
j_beg = nxi ;
Phi = zeros(Probl.problem_size,1) ;
for j = j_beg : n-1
Phi = Phi + PC.an(n-j+1)*fy(:,j+1) ;
end
Phi_n = St + PC.halpha2*(PC.a0(n+1)*fy(:,1) + zn_corr(:,n+1) + Phi) ;
yn0 = y_pred ; fn0 = f_pred ;
stop = 0 ; mu_it = 0 ;
while ~stop
yn1 = Phi_n + PC.halpha2*fn0 ;
mu_it = mu_it + 1 ;
if PC.mu == Inf
stop = norm(yn1-yn0,inf) < PC.mu_tol ;
if mu_it > 100 && ~stop
warning('MATLAB:fde12:NonConvegence',...
strcat('It has been requested to run corrector iterations until convergence but ', ...
'the process does not converge to the tolerance %e in 100 iteration'),PC.mu_tol) ;
stop = 1 ;
end
else
stop = mu_it == PC.mu ;
end
fn1 = f_vectorfield(t(n+1),yn1,Probl) ;
yn0 = yn1 ; fn0 = fn1 ;
end
y(:,n+1) = yn1 ;
fy(:,n+1) = fn1 ;
end
end

end

% =========================================================================
% =========================================================================
function z = FastConv(x, y)

r = length(x) ; problem_size = size(y,1) ;

z = zeros(problem_size,r) ;
X = fft(x,r) ;
for i = 1 : problem_size
Y = fft(y(i,:),r) ;
Z = X.*Y ;
z(i,:) = ifft(Z,r) ;
end

end

% =========================================================================
% =========================================================================
function f = f_vectorfield(t,y,Probl)

if isempty(Probl.param)
f = feval(Probl.fdefun,t,y) ;
else
f = feval(Probl.fdefun,t,y,Probl.param) ;
end

end

% =========================================================================
% =========================================================================
function ys = StartingTerm(t,ic)

ys = zeros(size(ic.y0,1),1) ;
for k = 1 : ic.m_alpha
ys = ys + (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;
end

end``````
LE_RF.m
``````function f=LE_RF(~,x)
%Output data must be a column vector
f=zeros(size(x));
%variables allocated to the variational equations
X= [x(4), x(7), x(10);
x(5), x(8), x(11);
x(6), x(9), x(12)];
%RF equations
f(1)=x(2)*(x(3)-1+x(1)*x(1))+0.1*x(1);
f(2)=x(1)*(3*x(3)+1-x(1)*x(1))+0.1*x(2);
f(3)=-2*x(3)*(0.98+x(1)*x(2));
%Jacobian matrix
J=[2*x(1)*x(2)+0.1, x(1)*x(1)+x(3)-1, x(2);
-3*x(1)*x(1)+3*x(3)+1,0.1,3*x(1);
-2*x(2)*x(3),-2*x(1)*x(3),-2*(x(1)*x(2)+0.98)];
%Righthand side of variational equations
f(4:12)=J*X;
end``````
LE_RF.m
``````function f=LE_RF(~,x)
%Output data must be a column vector
f=zeros(size(x));
%variables allocated to the variational equations
X= [x(4), x(7), x(10);
x(5), x(8), x(11);
x(6), x(9), x(12)];
%RF equations
f(1)=x(2)*(x(3)-1+x(1)*x(1))+0.1*x(1);
f(2)=x(1)*(3*x(3)+1-x(1)*x(1))+0.1*x(2);
f(3)=-2*x(3)*(0.98+x(1)*x(2));
%Jacobian matrix
J=[2*x(1)*x(2)+0.1, x(1)*x(1)+x(3)-1, x(2);
-3*x(1)*x(1)+3*x(3)+1,0.1,3*x(1);
-2*x(2)*x(3),-2*x(1)*x(3),-2*(x(1)*x(2)+0.98)];
%Righthand side of variational equations
f(4:12)=J*X;
end``````
F0_Lyapunov.m
``````function [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out)
ne=3;
ext_fcn=@LE_RF;
t_start=0;
h_norm=0.01;
t_end=10;
x_start=[0.1;0.1;0.1];
h=0.01;
q=0.75;
out=0;

% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Program to compute LEs as function of time (t) of autonomous systems of commensurate Fractional Order (FO)
% defined via Caputo's derivative;
%
% ***********************************************************************************************
% August 2022: Improved to overcome the problem of plot function of the last
% variants of Matlab
% ***********************************************************************************************
%
% author: Marius-F. Danca, 2018
% http://danca.rist.ro/
% email: danca@rist.ro
% The program uses a fast variant of the predictor-corrector Adams-Bashforth-Moulton,
% FDE12.m for FDEs, by Roberto Garrappa: https://goo.gl/XScYmi
% m-files required: FDE12, FO_Lyapunov.m and the function containing the extended system
% (see e.g. LE_RF.m).
% FO_Lyapunov.m, FDE12.m and LE_FO_RF.m must be in the same folder.
% FO_Lyapunov.m was developed, by modifying the program Lyapunov.m, by V.N. Govorukhin:
% https://goo.gl/wZVCtg
%
% How to use it:
% [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out);
%
% Input:
% ne - system dimension (eqs number);
% ext_fcn - function containing the extended system of FO;
% t_start, t_end - time span (FDE12);
% h - integration step;
% h_norm - step for Gram-Schmidt renormalization (Lyapunov algorithm), an important parameter:
%       experimentaly determined as an optimal choice if it is multiple of h, or equal;
% x_start - initial condition;
% q - 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.
%
% Example of use for the RF system:
% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Cite as:
%
% Marius-F. Danca and N. Kuznetsov, Matlab code for Lyapunov exponents of
% fractional order systems, International Journal of Bifurcation and Chaos,
% 28(05)(2018), 1850067
%
% Marius.-F. Danca, Matlab code for Lyapunov exponents of fractional-order systems,
% Part II: The non-commensurate case, International Journal of Bifurcation and Chaos,
% 31(12), 2150187, (2021)
%
hold on;
colors = 'grby';%can be expanded expand that to include any of the 8 color codes bcgkmrwy
% 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
x(1:ne)=x_start;
i=1;
while i<=ne
x((ne+1)*i)=1.0;
i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
LE=zeros(ne,1);
% Solutuion of extended ODE system of FO using FDE12 routine
[~,Y] = fde12(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
%       print results
if (mod(it,out)==0)
fprintf('%10.2f %10.4f %10.4f %10.4f\n',[t,LE']);
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;
%       display LEs
for i=1:ne
plot(t,LE(i),'.','Color',colors(i));
%plot(t,LE(i),'.','markersize',1,'Color',colors(i));%for thiny draw
end
end
xlabel('t','fontsize',16)
ylabel('LEs','fontsize',14)
set(gca,'fontsize',14)
box on;
line([0,t],[0,0],'color','k')
axis tight``````
F0_Lyapunov.m
``````function [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out)
ne=3;
ext_fcn=@LE_RF;
t_start=0;
h_norm=0.01;
t_end=10;
x_start=[0.1;0.1;0.1];
h=0.01;
q=0.75;
out=0;

% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Program to compute LEs as function of time (t) of autonomous systems of commensurate Fractional Order (FO)
% defined via Caputo's derivative;
%
% ***********************************************************************************************
% August 2022: Improved to overcome the problem of plot function of the last
% variants of Matlab
% ***********************************************************************************************
%
% author: Marius-F. Danca, 2018
% http://danca.rist.ro/
% email: danca@rist.ro
% The program uses a fast variant of the predictor-corrector Adams-Bashforth-Moulton,
% FDE12.m for FDEs, by Roberto Garrappa: https://goo.gl/XScYmi
% m-files required: FDE12, FO_Lyapunov.m and the function containing the extended system
% (see e.g. LE_RF.m).
% FO_Lyapunov.m, FDE12.m and LE_FO_RF.m must be in the same folder.
% FO_Lyapunov.m was developed, by modifying the program Lyapunov.m, by V.N. Govorukhin:
% https://goo.gl/wZVCtg
%
% How to use it:
% [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out);
%
% Input:
% ne - system dimension (eqs number);
% ext_fcn - function containing the extended system of FO;
% t_start, t_end - time span (FDE12);
% h - integration step;
% h_norm - step for Gram-Schmidt renormalization (Lyapunov algorithm), an important parameter:
%       experimentaly determined as an optimal choice if it is multiple of h, or equal;
% x_start - initial condition;
% q - 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.
%
% Example of use for the RF system:
% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Cite as:
%
% Marius-F. Danca and N. Kuznetsov, Matlab code for Lyapunov exponents of
% fractional order systems, International Journal of Bifurcation and Chaos,
% 28(05)(2018), 1850067
%
% Marius.-F. Danca, Matlab code for Lyapunov exponents of fractional-order systems,
% Part II: The non-commensurate case, International Journal of Bifurcation and Chaos,
% 31(12), 2150187, (2021)
%
hold on;
colors = 'grby';%can be expanded expand that to include any of the 8 color codes bcgkmrwy
% 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
x(1:ne)=x_start;
i=1;
while i<=ne
x((ne+1)*i)=1.0;
i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
LE=zeros(ne,1);
% Solutuion of extended ODE system of FO using FDE12 routine
[~,Y] = fde12(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
%       print results
if (mod(it,out)==0)
fprintf('%10.2f %10.4f %10.4f %10.4f\n',[t,LE']);
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;
%       display LEs
for i=1:ne
plot(t,LE(i),'.','Color',colors(i));
%plot(t,LE(i),'.','markersize',1,'Color',colors(i));%for thiny draw
end
end
xlabel('t','fontsize',16)
ylabel('LEs','fontsize',14)
set(gca,'fontsize',14)
box on;
line([0,t],[0,0],'color','k')
axis tight``````
...