Bifurcation Diagram for system of ODEs

We are trying to plot a bifurcation diagram for the Fitz-Hugh Nagumo model: This is our function - (these parameter values show chaotic behaviour)

We understand that we probably need to use loops to create the birfurcation diagram and were going to use the parameter a to vary against, but we can't figure out how to do this in MATLAB. Any help would be appreciated.

function dv=myode3(t,v)
dv=zeros(2,1);
c=0.757; 
b=0.124;
a=0.191;
f=0.1297;
l=(a/(2*pi*f))*cos(2*pi*f*t); 
e=0.12;
dv(1)=v(1)*((1-v(1))*(v(1)-b))-v(2)+l;
dv(2)=e*(v(1)-c*v(2));
end