Need Bifurcation Diagram from code please.

% Define the system of ODEs

function f = harmonic_oscillator(t, z, omega, gamma, gamma_drive)

x = z(1);

x_dot = z(2);

f = [x_dot; -2 * gamma * x_dot - omega^2 * x + gamma_drive * cos(1.25 * t)];

end

% Set the parameters

omega = 1.25;

gamma_values = [0.2, 0.3, 0.315, 0.37, 0.5, 0.8];

t_span = [0, 400];

% Loop through different gamma values

for gamma = gamma_values

% Set the initial conditions

initial_conditions = [1; 0];

% Define the function for the ODE solver

ode_func = @(t, z) harmonic_oscillator(t, z, omega, gamma, gamma);

% Use the ODE solver to obtain the solution

[t, sol] = ode45(ode_func, t_span, initial_conditions);

% Extract the solution

x = sol(:, 1);

x_dot = sol(:, 2);

% Plot the numerical solution

figure('Position', [100, 100, 1200, 400]);

subplot(1, 2, 1);

plot(t, x, 'DisplayName', 'x(t)');

xlabel('Time');

ylabel('Displacement');

title(['Numerical Solution for Gamma = ', num2str(gamma)]);

legend;

% Plot the phase portrait

subplot(1, 2, 2);

plot(x, x_dot, 'DisplayName', 'Phase Portrait');

xlabel('x');

ylabel('x_dot');

title(['Phase Portrait for Gamma = ', num2str(gamma)]);

legend;

drawnow;

end

This is my code and I am not sure on how to do the bifurcation plots now.