How to select one period in non-linear system response?

Hi all,

I have a non-linear SDOF system excited by a sine wave with increasing frequency:

y = y0 *sin(2*pi*fe*t)

and I'm using ode15s to obtain the displacement of the mass. The result is a sine wave with a small non-linear transient response. I need to select one piece of this response (which is ideally in steady state), namely fourth period for each response. Below you can see a simplified mock-up of my code. I'm not sure the method I'm using to select the cycle is correct in the long haul (it seems to work only for excitation frequencies below like 50).

clc
clear all
close all


% Linear Parameters________________________________________________________
K=5000;% linear stiffness
M=1; % system's mass

%Input Parameters__________________________________________________________
y0 = 1; % amplification of base input


% Non-linear Parameters__________________________________________________________
Fs=0.5; % mu_static
s0=2000;% sigma 
Fn=400; %normal force
mu = 0.3;


% Excitation sweepsinus____________________________________________________
fe = 1:0.1:10;          %excitation 
n_freq = length(fe);    %number of simulation iterations

XR = [];
z0=zeros(3,1); % initialisation
cycle_i = 4;    %index of period to sample
cycle_f = cycle_i + 1;

for kk=1:n_freq

    fs = 64*fe(kk);         %sampling frequency
    dt = 1/fs;              %time step
    np = 20;                %number of periods
    T = 1/fe(kk);           %period
    T_s = T*np;             %simulation time
    N = T_s*fs +1 ;        %samples for each simulation time
    t= (0:N-1)'*dt;         %time vector of simulation time
    df=1/T_s;               %frequency step
    f = 0:df:fs;            %frequency vector
    
%_________________________SOLVING ODE FOR fe_______________________________
    

    y = y0*sin(2*pi*fe(kk)*t);     %excitation input
    y_dot = 2*pi*fe(kk)*y0*cos(2*pi*fe(kk)*t);
    
    [t,z] = ode15s('system',(0:dt:T_s),z0,[],s0,Fn,K,M, mu, fe(kk), y0 );

    x=z(:,1);
    v=z(:,2);
    
    % Period sampling for classical continuation
    initial_sample = cycle_i * T;   %time at which the fourth period starts
    final_sample = cycle_f * T;   %time at which the fifth period starts i.e fourth period ends
    
    idx_i = find(abs(t-initial_sample)<0.001);  %index for which fourth period starts
    idx_f = find(abs(t-final_sample)<0.001);    %index for which fourth period ends
    
    z0 = z(idx_f,:);                            %new initial conditions 
    
    sample = x(idx_i:idx_f,:);                  %sample between 4th and 5th period
    sample_time = t(idx_i:idx_f,:);             %time of new sample
 
    figure(1)
    subplot(2,1,1)
    plot(t, x, t(idx_i),sample(1,:),'o', t(idx_f),sample(end,:),'o')
    subplot(2,1,2)
    plot(sample_time, sample)
    
end

So in the end, I'm interested in the signal between 4*T and 5*T for each excitation frequency fe (the plots inside the loop are to monitor the sampling). Is this a good way to do this, or am I overcomplicating the matters? Also, how would you check if the sampled part is in steady state?

Thank you in advance!