I want to plot the bifurcation diagram of forced Duffing oscillations versus excitation frequency such as for Logistic Map ?

code for Logistic Map :

clear; clc;

% Parameters
N = 1000; % Number of iterations
mu_range = linspace(0, 4, 1000); % Range of mu values

% Preallocate arrays
x = zeros(N, length(mu_range));
lyapunov = zeros(1, length(mu_range));

% Loop through mu values
for i = 1:length(mu_range)
    mu = mu_range(i);
    x(1,i) = 0.5; % Initial condition
    
    % Iterate Logistic map
    for j = 2:N
        x(j,i) = mu * x(j-1,i) * (1 - x(j-1,i));
    end
    
    % Calculate Lyapunov exponent
    lyapunov_sum = 0;
    for j = 2:N
        lyapunov_sum = lyapunov_sum + log(abs(mu - 2 * mu * x(j,i)));
    end
    lyapunov(i) = lyapunov_sum / (N - 1);
end

% Plot bifurcation diagram
figure;
plot(mu_range, x(end-200:end,:), '.', 'MarkerSize', 1);
xlabel('\mu');
ylabel('x');
title('Bifurcation Diagram of Logistic Map');

% Plot Lyapunov exponent
figure;
plot(mu_range, lyapunov);
xlabel('\mu');
ylabel('\lambda');
title('Lyapunov Exponent vs \mu');

code for forced duffing oscillations :

clear; clc;

% Parameters
N = 1000; % Number of iterations
Omega_range = linspace(0.5, 2.5, 1000); % Range of Omega values

% Preallocate arrays
x = zeros(N, length(Omega_range));
lyapunov = zeros(1, length(Omega_range));

% Loop through Omega values
for i = 1:length(Omega_range)
    Omega = Omega_range(i);
    F = 1.0; % For simplicity, assume a constant force
    
    % Solve the forced Duffing equation using a numerical solver (e.g., ode45)
    [~, sol] = ode45(@(t, x) duffing(t, x, F, Omega), linspace(0, 100, N), [0.1; 0]);
    x(:, i) = sol(:, 1); % Store position
    
    % Calculate Lyapunov exponent
    lyapunov_sum = 0;
    for j = 2:N
        lyapunov_sum = lyapunov_sum + log(norm(duffing_jacobian(sol(j,:))));
    end
    lyapunov(i) = lyapunov_sum / (N - 1);
end

% Plot bifurcation diagram
figure;
plot(Omega_range, x(end-200:end,:), '.', 'MarkerSize', 1);
xlabel('\Omega');
ylabel('x');
title('Bifurcation Diagram of Forced Duffing Oscillation');

% Plot Lyapunov exponent
figure;
plot(Omega_range, lyapunov);
xlabel('\Omega');
ylabel('\lambda');
title('Lyapunov Exponent vs \Omega');

% Duffing equation
function dx = duffing(t, x, F, Omega)
    gamma = 1;
    delta = 0.4;
    omega = 0.5;
    a = 0.25;
    dx = [x(2); gamma*cos(omega*t) - (delta*x(2) + x(1)*(x(1)^2 - 1)*(x(1)^2 - a) - F*x(1)*cos(Omega*t))];
end

% Jacobian of Duffing equation
function J = duffing_jacobian(x)
    gamma = 1;
    delta = 0.4;
    a = 0.25;
    Omega = 1.5;
    J = [0, 1; -(delta + 6*x(1)^5 - 6*a*x(1)^3 - x(1)*Omega*sin(Omega)), -delta];
end

Could you please provide more details on what specifically is not working correctly with the code for Duffing oscillations? This will help me better understand the issue and provide an appropriate solution.