ClearAll["Global`*"]​​​​M=1;​​K=1;​​C=0.1;​​C1=0.1;​​Q=1;​​P=1;​​R0=1;​​Xp=1;​​Uinf=1;​​b0=1;​​​​eta1=η[τ][[1]];​​eta2=η[τ][[2]];​​​​eqs={Mη''[τ]+Kη[τ]+(Cη'[τ]-C1z'[τ])==Q(η[[2]]+η'[[1]]/Uinf+b0η'[[2]])^3+PV[τ],V'[τ]+R0V[τ]+Xp(η'[[1]]-z'[τ])==0};​​ClearAll[l1,l2]​​l1={};​​l2={};​​Table[NDSolve[{eqs}~Join~{η[0]=={0,0},η'[0]=={0,0},z[0]==0,z'[0]==0,V[0]==0},WhenEvent[η'[[1]][τ]==0,If[τ>100,{AppendTo[l1,cp],AppendTo[l2,η[τ]]}]]],{η[τ],z[τ],V[τ]},{τ,0,120}],{cp,0.1,2,0.01}];​​​​data=Transpose[{l1,l2}];​​ListPlot[data,Frame->True,FrameStyle->{Black,Medium},Axes->False,FrameLabel->{"β1","η"},PlotStyle->{PointSize[0.001],Opacity[0.7],Red}]