How might i plot bifurcation diagram

I tried to plot a bifurcation diagram about the hindmarsh-rose model. I made its nonlinear-equations compute with fsolve in terms of the parameter u in the equations. There is no here problem.

But I couldn't plot the bifurcation diagram.

Can anyone help me about this case? Thanks for your help

My Code is here :

 function F = Hindmarsh_Rose_84_Equation(x,u)
X=x(1);
Y=x(2);
Z=x(3);
F = [Y-X.^3+3*X.^2-Z+u;
1-5*X.^2-Y;
0.001*(4*(X-(-1.618))-Z);
];
end
hindmarsh_rose = @Hindmarsh_Rose_84_Equation;
u =0:0.01:5;
x0 = [0.1,0.1,0.1];
for i=1:length(u)
%    options = optimoptions('fsolve','Display','iter'); 
x =fsolve(hindmarsh_rose,x0,[],u(i));

    disp('u');
    disp(u(i));
    disp('x,y,z');
    disp(x(1:3));
   end

   plot(u(i),x(2),'k.','MarkerSize',2);

end