Nonlinear frequency response – Mathematica Stack Exchange

I am trying to reproduce the result of this paper, namely Figure 10 which depicts the solution of equation 64 given as,

enter image description here

qmax is the steady state solution of the oscillatory equation given. I tried to reproduce the result using the code below:

data = {Vac -> 0.01 Cos((CapitalOmega) t), Vdc -> 100, c -> 0};
 eq = q''(t) + c*q'(t) + 864.82*q(t) + 21.78 q(t)^2 + 142.58 q(t)^3     == 
  45.31 ((2 Vac)/Vdc + Vac^2/Vdc^2) + (1 + (2 Vac)/Vdc + Vac^2/
  Vdc^2)*(123.51 q(t) + 74.09 q(t)^2 + 1424.72 q(t)^3 - 
  122.45 q(t)^4);
sol = ParametricNDSolve({eq /. data, q(0) == 0, q'(0) == 0}, 
 q, {t, .1, 100}, {(CapitalOmega)});
lst={#, Max(q(#)(Range(80, 100, .1)) /. sol)} & /@ Range(20, 30, .1);

ListPlot(lst);

As you can see the figure is different. I appreciate if you can help in getting the same figure (figure 10), shown below

enter image description here

Regards,