How I cant Plot the Spectrum of Lyapunov exponents of system versus the parameter c∈ [0, 200], with a = 10, b = 28 ?

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