Greetings, all. I trust you're doing well. My intention is to generate a spectrum depicting the Lyapunov exponents of the system concerning the parameter c within the range of \( [0, 200] \), given ( a = 10 ) and ( b = 28 ). Fig(11.png).The system equations are as follows:

x'= -a * x + b * y - y * z;

y' = x + z * x;

z' = -c * z + y^2;

I'm employing specific codes to calculate the Lyapunov exponents for c = 2 Fig(22.png).

codes :

%%%%%%%%% system_equations %%%%%%%%%%%%%%%%%%%

% Define the system equations

function f = system_equations(t, X)

% Define the parameters

a = 10;

b = 28;

c = 2;

% Extract variables

x = X(1);

y = X(2);

z = X(3);

% Extract the Jacobian matrix Y from the extended state vector X

Y = [X(4), X(7), X(10);

X(5), X(8), X(11);

X(6), X(9), X(12)];

% Initialize the output vector

f = zeros(9, 1);

% System equations

f(1) = -a * x + b * y - y * z;

f(2) = x + z * x;

f(3) = -c * z + y^2;

% Jacobian matrix

Jac = [-a, b - z, -y;

1 + z, 0, x;

0, 2 * y, -c];

% Variational equation

f(4:12) = Jac*Y;

end

%%%%%%%%% systemequations_Lyapunov %%%%%%%%%%%%%%%%%%%

function [Texp, Lexp] = systemequations_Lyapunov(n, rhs_ext_fcn, fcn_integrator, tstart, stept, tend, ystart, ioutp)

clc;

close all;

n = 3;

rhs_ext_fcn = @system_equations;

fcn_integrator = @ode45;

tstart = 0;

stept = 0.1;

tend = 1000;

ystart = [7.2 7.8 2.3];

ioutp = 10;

n1 = n;

n2 = n1*(n1 + 1);

% Number of steps.

nits = round((tend - tstart)/stept);

% Memory allocation.

y = zeros(n2, 1);

cum = zeros(n1, 1);

y0 = y;

gsc = cum;

znorm = cum;

% Initial values.

y(1:n) = ystart(:);

for i = 1:n1

y((n1 + 1)*i) = 1.0;

end

t = tstart;

% Main loop.

for iterlya = 1:nits

% Solutuion of extended ODE system.

[T, Y] = feval(fcn_integrator, rhs_ext_fcn, [t t + stept], y);

t = t + stept;

y = Y(size(Y, 1),:);

for i = 1:n1

for j = 1:n1

y0(n1*i + j) = y(n1*j + i);

end

end

% Construct new orthonormal basis by Gram-Schmidt.

znorm(1) = 0.0;

for j = 1:n1

znorm(1) = znorm(1) + y0(n1*j+1)^2;

end

znorm(1) = sqrt(znorm(1));

for j = 1:n1

y0(n1*j + 1) = y0(n1*j + 1)/znorm(1);

end

for j = 2:n1

for k = 1:(j - 1)

gsc(k) = 0.0;

for l = 1:n1

gsc(k) = gsc(k) + y0(n1*l + j)*y0(n1*l + k);

end

end

for k = 1:n1

for l = 1:(j - 1)

y0(n1*k + j) = y0(n1*k + j) - gsc(l)*y0(n1*k + l);

end

end

znorm(j) = 0.0;

for k = 1:n1

znorm(j) = znorm(j) + y0(n1*k+j)^2;

end

znorm(j)=sqrt(znorm(j));

for k = 1:n1

y0(n1*k + j) = y0(n1*k + j)/znorm(j);

end

end

% Update running vector magnitudes.

for k = 1:n1

cum(k) = cum(k) + log(znorm(k));

end

% Normalize exponent.

for k = 1:n1

lp(k) = cum(k)/(t - tstart);

end

% Output modification.

if iterlya == 1

Lexp = lp;

Texp = t;

else

Lexp = [Lexp; lp];

Texp = [Texp; t];

end

for i = 1:n1

for j = 1:n1

y(n1*j + i) = y0(n1*i + j);

end

end

end

% Show the Lyapunov exponent values on the graph.

str1 = num2str(Lexp(nits, 1));

str2 = num2str(Lexp(nits, 2));

str3 = num2str(Lexp(nits, 3));

pl = plot(Texp, Lexp,'LineWidth',2);

pl(1).Color = 'b';

pl(2).Color = 'g';

pl(3).Color = 'r';

legend('\lambda_1', '\lambda_2','\lambda_3')

legend('Location','northeast')

title('Dynamics of Lyapunov Exponents');

text(205, 1.5, '\lambda_1 = ','Fontsize',20);

text(220, 1.68, str1,'Fontsize',20);

text(205, -1, '\lambda_2 = ','Fontsize',20);

text(220, -0.73, str2,'Fontsize',20);

text(205, -13.8, '\lambda_3 = ','Fontsize',20);

text(220, -13.5, str3,'Fontsize',20);

xlabel('Time');

ylabel('Lyapunov Exponents');

axis([-1,300, -16,6]);

set(gca,'FontSize',20)

end

% End of plot