How to get the frequency response with the unstable branches for a nonlinear driven system?

Here's a partial answer, because I suspect the OP's code doesn't correspond perfectly with their figure, resulting in a more complicated situation.

The main idea is to make a stroboscopic map that advances the system by one period. This can be treated as a discrete-time dynamical system. Equilibria (corresponding to limit cycles of the original system) can be found by FindRoot, including unstable cycles. Finally we can implement a crude continuation technique to track an equilibrium across the parameter, by using the previous answer as the initial guess for FindRoot.

First, a helper to find extrema:

FindExtrema[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}] :=
  Reap[NDSolve[{y'[t] == Evaluate[D[if[t], t]],
  WhenEvent[y'[t] == 0, Sow[{t, y[t]}]],
  y[tmin] == if[tmin]}, y[t], {t, tmin, tmax}]][[2]] /. {x_List} :> x;

FindExtrema[if_InterpolatingFunction] := FindExtrema[if, if["Domain"][[1]]];

The model is set up as in the original post (not copied here), except I changed the second equation to use y[t] instead of x'[t] to avoid the need for Method -> {"EquationSimplification" -> "Residual"}]:

eqns = {
  x'[t] == y[t],
  y'[t] == -(1/2 y[t]^2 (3 - δ/x[t]^3 (1 + δ/x[t]^3)^(-4/3) - 
    3 (1 + δ/x[t]^3)^(-1/3)) + intf[x[t], z[t]] - deltap - p0*Sin[ω*t])/x[t]/(1 - (1 + δ/x[t]^3)^(-1/3)), 
  z'[t] == z[t]*(μ (x[t]^2*z[t]^-2 - x[t]^-4*z[t]^4))/(3 ηb*(1 - (x[t]^-4*z[t]^4 + 2 x[t]^2*z[t]^-2 - 3)/β))};

Now, define the stroboscopic map F:

F[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] := (
  sol = NDSolve[{eqns, {x[0] == x0, y[0] == y0, z[0] == z0}},
    {x, y, z}, {t, 0, 2 π/ω}][[1]];
 {x[2 π/ω], y[2 π/ω], z[2 π/ω]} /. sol
);

Next a function to find fixed points of the map:

findEq[{x0i_?NumericQ, y0i_?NumericQ, z0i_?NumericQ}, opts___] :=
  {x0, y0, z0} /. 
  FindRoot[F[{x0, y0, z0}] == {x0, y0, z0}, {{x0, x0i}, {y0, y0i}, {z0, z0i}}, opts]

In action:

ω = 2.5;
eq = findEq[{1.2, 0.1, 1.2}]
(* {1.23676, 0.169015, 1.20327} *)
Plot[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 2 π/ω}]

Since FindRoot may throw some FindRoot::lstol warnings, we can check whether the equilibrium we found is good:

F[eq] - eq
(* {-8.50313*10^-9, -4.05657*10^-8, 4.94302*10^-10} *)

Close enough!

We can also check the stability of a cycle with a finite-difference approximation to the Jacobian matrix of F:

j := {(F[eq + {ϵ, 0, 0}] - F[eq - {ϵ, 0, 0}])/(2 ϵ),
      (F[eq + {0, ϵ, 0}] - F[eq - {0, ϵ, 0}])/(2 ϵ),
      (F[eq + {0, 0, ϵ}] - F[eq - {0, 0, ϵ}])/(2 ϵ)};

ϵ = 10^-5;
Chop[Eigenvalues[j]]
(* {0.852387, 0.0131683 + 0.598595 I, 0.0131683 - 0.598595 I} *)

All real parts less than one in magnitude indicates stability. ϵ needs to be small enough to make this a good linearization of F' but not too small, to avoid numerical roundoff errors.

Now we'll make four tracks of the equilibria with different starting values of ω:

ω = 2.5;
ics = findEq[{1.2367572037134027`, 0.169014899645155`, 1.2032719485392172`}];
Clear[ω];

res1 = Table[
   ics = eq1[ω] = findEq[ics];
   {ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
   , {ω, 2.5, 3.27, 0.01}];

ω = 3.5;
ics = findEq[{1.7, -3, 1.9}];
Clear[ω];

res2 = Table[
   ics = eq2[ω] = findEq[ics];
   {ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
   , {ω, 3.5, 2.7, -0.01}];

ω = 3.0;
ics = findEq[{0.81, -1.31, 1.51}, DampingFactor -> 0.1];
Clear[ω];

res3 = Table[
   ics = eq3[ω] = findEq[ics];
   {ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
   , {ω, 3.0, 3.27, 0.01}];

ω = 3.0;
ics = findEq[{0.81, -1.31, 1.51}, DampingFactor -> 0.1];
Clear[ω];

res4 = Table[
   ics = eq4[ω] = findEq[ics];
   {ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
   , {ω, 3.0, 2.5, -0.01}];

res3 and res4 are an unstable cycle. Getting the initial point for these tracks is the hardest part and involves some trial and error.

Plotting the results:

ListPlot[{res1, res2, res3, res4}]

Note the discrepancy between these results and OP's figure. The unstable branch doesn't connect up with the upper stable branch, so I suspect there is a lot more going on in this system. But without a reference to the source of the figure, it's hard to know how to proceed.

By the way, here's an alternative equilibrium tracking code that uses linear extrapolation of the last two answers to get an even better initial guess for FindRoot:

ω = 2.5;
ics′ = ics = findEq[{1.2367572037134027`, 0.169014899645155`, 1.2032719485392172`}];
Clear[ω];

res1′ = Table[
   eq1[ω] = findEq[2 ics - ics′];
   ics′ = ics; ics = eq1[ω];
   {ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
   , {ω, 2.5, 3.27, 0.01}];