Bifurcation analysis of a ODE system

clear;
% Parameters
bifurcationParam = linspace(0, 10, 100);  % Range of parameter values for the bifurcation analysis
a = rand;    % Fixed system parameter a
c = rand;    % Fixed system parameter c
d = rand;    % Fixed system parameter d
% Define the ODE system
dxdt = @(t, x, y, z) a * (y - x);
dydt = @(t, x, y, z, bifParam) x - bifParam * x * z + c * y;
dzdt = @(t, x, y, z) x * y + x^2 * z - d * z;
% Bifurcation analysis
bifurcationData = zeros(length(bifurcationParam), 100);
x0 = [1; 1; 1];
tspan = [0 100];
for i = 1:length(bifurcationParam)
    currentParam = bifurcationParam(i);
    odeSystem_i = @(t, X) [
        dxdt(t, X(1), X(2), X(3));
        dydt(t, X(1), X(2), X(3), currentParam);
        dzdt(t, X(1), X(2), X(3));
        ];
    
    % Compute solution only at specific time points
    [~, y] = ode113(odeSystem_i, tspan, x0);
    sampleIndices = round(linspace(1, length(y), 100));
    bifurcationData(i, :) = y(sampleIndices, 1)';
end
% Plot bifurcation diagram
figure;
plot(bifurcationParam, bifurcationData, 'r.', 'MarkerSize', 10);
xlabel('Parameter y');
ylabel('x');
title('Bifurcation Diagram');

Here b I have chosen to be a bifurcation parameter. But it does not give proper bifurcation. Can anyone help.

[YOUSSEF El MOUSSATI]

It seems like there might be an issue with how the bifurcation parameter b is used in the system dynamics. In the code provided, b is not defined or used in the ODE system. Instead, bifParam is used as the bifurcation parameter in the dydt function.

To fix this issue, you can replace b with bifParam in the dydt function definition, like this:

    clear;

% Parameters bifurcationParam = linspace(0, 10, 100); % Range of parameter values for the bifurcation analysis a = rand; % Fixed system parameter a c = rand; % Fixed system parameter c d = rand; % Fixed system parameter d

% Define the ODE system dxdt = @(t, x, y, z) a * (y - x); dydt = @(t, x, y, z, bifParam) x - bifParam * x * z + c * y; dzdt = @(t, x, y, z) x * y + x^2 * z - d * z;

% Bifurcation analysis bifurcationData = zeros(length(bifurcationParam), 100); x0 = [1; 1; 1]; tspan = [0 100];

for i = 1:length(bifurcationParam) currentParam = bifurcationParam(i); odeSystem_i = @(t, X) [ dxdt(t, X(1), X(2), X(3)); dydt(t, X(1), X(2), X(3), currentParam); dzdt(t, X(1), X(2), X(3)); ];

    % Compute solution only at specific time points
    [~, y] = ode113(odeSystem_i, tspan, x0);
    sampleIndices = round(linspace(1, length(y), 100));
    bifurcationData(i, :) = y(sampleIndices, 1)';
end

% Plot bifurcation diagram figure; plot(bifurcationParam, bifurcationData, 'r.', 'MarkerSize', 10); xlabel('Parameter b'); ylabel('x'); title('Bifurcation Diagram');