Sweep-up for Duffing equation

Question

I would like to calculate a sweep-up for the Duffing equation. To be more precise, I'm looking for a solution to the ode

$\qquad x''(t)+0.1x'(t)+x(t) +0.01 x^3(t) = 0.8 \cos(\eta t) $

Let's say for the initial conditions $x(0)=x'(0)=0$ and $\eta \in [0.8, 1.2]$ with increment size $0.001$ and $t \in [0,2\cdot 100 \frac{\pi}{\eta}]$. But instead of using the same initial conditions for every $\eta$, I would like to choose new initial conditions $x_{j+1}(0)=x_j(2\cdot 100 \frac{\pi}{\eta})$ and $x_{j+1}'(0)=x_j'(2\cdot 100 \frac{\pi}{\eta})$ every time $\eta$ increases, where $x_j$ is the solution of the previous computation with the previous $\eta$. As output I would like to have the maximum value of $x$ for every $\eta$ neglecting the transients.

I tried it using Table, which yields reasonable results but takes some time to calculate. Here is the code

numbper = 100; (*number of periods*)
sol = ParametricNDSolveValue[{x''[t] + 1/10 x'[t] + x[t] + 1/100 x[t]^3 == 
8/10 Cos[η t], x[0] == x0, x'[0] == v0}, {x, x'}, {t, 0, numbper*2 Pi/η},
{x0, v0, η}];

xstart = sol[0, 0, 0.79][[1]][numbper*2 Pi/0.79];
vstart = sol[0, 0, 0.79][[2]][numbper*2 Pi/0.79];
dη = 0.001;
sweepup = Table[{
 NMaxValue[{sol[xstart, vstart, η][[1]][
    t], (numbper - 1)*2 Pi/η <= t <= numbper*2 Pi/η}, 
  t], xstart = 
  sol[xstart, vstart, η][[1]][(numbper - 1)*2 Pi/η];
 vstart = sol[xstart, vstart, η][[2]][numbper*2 Pi/η];}
, {η, 0.8, 1.2, dη}];

I'm now wondering what I could do to speed things up.

Answer 1

FoldList won't give a big raise to the efficiency (the most time consuming part is equation solving), but it does make the code conciser:

sol[eta_] := 
  ParametricNDSolveValue[{x''[t] + 1/10 x'[t] + x[t] + 1/100 x[t]^3 == 8/10 Cos[eta t], 
    x[0] == x0, x'[0] == v0, 
    WhenEvent[x'[t] < 0, Sow[x[t], eta], "DetectionMethod" -> "Interpolation"]}, 
     {x[100], x'[100]}, {t, 0, 100}, {x0, v0}];
etalst = Range[8/10, 1, 1/100];
{iclst, maxlst} = Reap@FoldList[sol@#2 @@ # &, {0, 0}, etalst];

maxValue = Max /@ ({Partition[iclst\[Transpose] // First, 2, 1], maxlst}\[Transpose])

(* {2.77457, 3.50534, 2.82074, 2.40488, 2.96662, 3.81174, 4.24709, 4.08059, 3.62443, 
2.8937, 3.46553, 4.14776, 4.67875, 4.97169, 4.68464, 3.92852, 4.15529, 4.71447, 5.13502, 
5.53676, 5.53677} *)