Bifurcation diagram plot command?

I am trying to produce a bifurcation diagram. I am new to the bifurcation world so i need a little help I have been able to produce the file that runs, however i cant plot the graph.

    for ilop=20000:100:30000
    global  ms Ixx ks1 ks2 coff1 coff2 lf lr pha w Fo pc lp tf

        close all
    clear
    %clc   %CLC clears the command window and homes the cursor

    load stiffness

    %System parameters .
    Fo=0.1;   w=6.9;

    ms=730; Ixx=2460; ks1=17500; ks2=17500; 
    coff1=784; coff2=784; 
    lf=1.011; lr=1.803;
    pha=0.3*pi;

    ks1=ks1+100;

    pc=0;  lp=1;
    xs0=[0  0  0  0  0  0];   stsz=0.001;

    i0=0;i1=0;i2=0;i3=0;i4=0;i5=0;i6=0;ig=0;   

    n=6000; tf=((2*pi/w)*n);            st=2000;

    Istep=2*pi/(w*st);  lep=0:Istep:(tf-Istep); lp1=length(0: (2*pi/(w*st)):(2*pi/w)-(2*pi/(w*st)));  

    lop=round(length(lep)/lp1);  

    xs0=[0.0  0.0  0.0  0.0  0.0  0.0]; 

    %Solver parameters setting

    options=odeset('AbsTol',1*exp(-12),'RelTol',1*exp(-12),'InitialStep',0.0000000001,'MaxStep',0.001);

    [t x]=ode45(@NLQCM1,[0:(2*stsz):tf],[xs0(1) xs0(2) xs0(3) xs0(4) xs0(5) xs0(6) ]);
    %ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %round(Nperiod); 
    rtf=round(3.5*length(t)/5);
    xT=[t x];
    Xp=xT(rtf:length(t),: );
    itl=0;

    for tL=  1   : lp1   :   length(Xp)-lp1
        itl=itl+1;
        se=[ks1 Xp(tL,:)];    
        xpm(itl,:)=se;
    end
    %Xt=[t x];
    load HCPM
    Xtot=[Xtot; xpm(:,1:6)];

    save HCPM Xtot ;
    save Stiffness ks1 ks2;
    end

The function file is

function [xs F]=HM(t,x)

global i ms Ixx ks1 ks1n ks2 ks2n coff1 coff2  lf lr  w Fo 

%state space vector
x1=x(1); x2=x(2); x3=x(3); x4=x(4);  x5=x(5); x6=x(6); 

i+1;

xs=[x1 x2  x3  x4  x5 x6]; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%State space model equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     u=Fo*sin(w*t); 

       Fo=0.05;   w=6.9;

ms=730; Ixx=2460; ks1=17500; ks2=17500; ks1n=34000; ks2n=34000;

coff1=784; coff2=784; 

lf=1.011; lr=1.803;pha=0.3*pi;  

      xg1=Fo*sin(w*t); xg2=Fo*sin(w*t-pha); xg1d=Fo*w*cos(w*t); xg2d=Fo*w*cos(w*t-pha);

      Xr=x(1)+lr*x(3);
      Xf=x(1)-lf*x(3);
      Xrd=x(2)+x(4)*lr;
      Xfd=x(2)-x(4)*lf;

       %if(i>=1)
       %    i=i+1
       %F(i)=u ;
      % end

       xs(1)=x2;

       xs(2)=(-(coff2*Xrd+coff1*Xfd+(ks2*Xr+ks2n*Xr^3)+(ks1*Xf+ks1n*Xf^3))/ms)+((coff2*xg2d+coff1*xg1d+(ks2*xg2+ks2n*xg2^3)+(ks1*xg1+ks1n*xg1^3))/ms);

       xs(3)=x4;

       xs(4)=(-(coff2*Xrd*lr-coff1*Xfd*lf+(ks2*Xr+ks2n*Xr^3)*lr+(-ks1*Xf-ks1n*Xf^3)*lf)/Ixx)+((coff2*xg2d*lr-coff1*xg1d*lf+(ks2*xg2+ks2n*xg2^3)*lr+(-ks1*xg1-ks1n*xg1^3)*lf)/Ixx);

       xs(5)=Fo*cos(w*t); 

       xs(6)=Fo*sin(w*t);
      xs=xs';

   end