ODE Bifurcation Diagram for loop

Hello everyone. I am having a lot of trouble with producing a bifurcation diagram for my project.

I have this set of equations :

dy(1) = y(2);

dy(2) = -y(1) -y(2)/Q - f0*sin(y(1) + w*t) + a;

I have to plot y(1) vs f0. I should pick a few values-points of y(1) that are equally spaced by 2*pi/w and plot them against the f0 (the value that the function used when it solved the ode). Below I have the code I've written so far and below the code I have a solution that I can't apply to my code. :

Solver :

clear all;

QWF; % call variables, see below

initial=[0 0];

time= 0:2*pi/w*10^-2:1000*2*pi/w;

[T,Y] = ode45(@myfunfinal,time,initial);

l = length(Y(:,1));

for n = 1:5

   hold all;

    plot(T(l - 100*n),Y(l-100*n,1),'x');

end

hold off;

The function that is meant to be solved is this

function dy = myfunfinal(t,y)

dy = zeros(2,1);

QWF;

a = - w/Q;

dy(1) = y(2);

dy(2) = -y(1) -y(2)/Q - f0*sin(y(1) + w*t) + a;

Where QWF is the .m file that has the constants in it.

Q = 20; w = 1; f0 = 0.1;

So what I need to do is to find a way to plot those values of y(1) vs many values of f0.

For this reason a for loop that will change the value of f0 in the function, then have the run the ode code and plot the results simultaneously is essential, however I can't get my head around it.

What I've thought so far is this.

for n = 1:990 f0 = 0.1 + 0.01*n end

or f0 = [0.1:0.01:10] for n = 1:990 f0(n) end