Bifurcation diagram of a 5-dimensional differential model

Hi all, I need to a plot a bifurcation diagram of a 5-dimensional differentioan model given below with respect to one of the parameters as a 'bifurcation parameter',(say 0.01< gamma<0.5) but I can't find an easy way to plot it. I am posting this question for the second time, any one can help me? Thanks

% % Solving a model using ode45
function PMWZF 
global r1 r2 r3 K1 K2 K3 eta1 eta2 eta3 eta4 eta5 eta6 alpha1 alpha2 P0 M0 W0 Z0 F0 a1 a2 b1 b2 alpha beta gamma delta dz df
% % Parameters 
r1=1;r2=0.3;r3=0.3;K1=20;K2=10;K3=10;eta1=0.01;eta2=0.01;eta3=0.1;eta4=0.1;
eta5=0.01;eta6=0.01;alpha1=0.11;alpha2=0.1;a1=0.1;a2=0.1;b1=0.01;b2=0.01;
alpha=0.024;beta=0.03;gamma=0.04;delta=0.04;dz=0.1;df=0.05;
% % initial conditions
P0 = 15;
M0 = 5;
W0=  8;
Z0 = 10;
F0 = 5;
init= [P0; M0; W0;Z0; F0];
tspan =[0 1600];
options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3 4 5]);
function dv = fun_PMWZF(t,v)
dv = zeros(5,1);
P = v(1); 
M = v(2); 
W = v(3); 
Z = v(4);
F = v(5); 
dv(1) = r1*P*(1- P/K1)- eta1*P*M-eta2*P*W- alpha1*P*Z- a1*P*F/(1+alpha*P+beta*M+gamma*W+delta*Z);
dv(2) = r2*M*(1- M/K2)- eta3*P*M-eta4*M*W;
dv(3)=  r3*W*(1-W/K3)- eta5*P*W-eta6*M*W;
dv(4) = alpha2*P*Z-a2*Z*F/(1+alpha*P+beta*M+gamma*W+delta*Z) - dz*Z ;
dv(5) = (b1*P*F+b2*Z*F)/(1+alpha*P+beta*M+gamma*W+delta*Z) -df*F;
end
[t,v] = ode45(@fun_PMWZF,tspan,init,options)
set(groot, 'DefaultAxesFontSize', 9)
plot(t,v,'LineWidth', 2)
xlabel('Time')
ylabel('Population')
title('PMWZF MODEL')
legend ('P','M','W','Z','F')
end