Frequency sweep for DAE question

%%
close all;
clear all;
clc;
format long;

k1=25000;
knl_1=300;
knl_3=2000000;
m=0.1;
c=1;

w0=20;
kAmp = 0.5;

M=[0,0,0; 0,1,0; 0,0,1];
opt = odeset('mass', M, 'RelTol',1.e-14);

a = w0-5;
b = w0+5;
np = 200;
f = linspace(a,b,np);
omega_rng = f;
y0 = [0; 0; 0]; % initial condition

for mm=1:length(omega_rng)    
    format long;
    omega=omega_rng(mm);
    w = omega;    
    period=2*pi/omega;
    step=period/100;
    endperiod=1000*period;    
    tspan=step: step: endperiod;    
    [t,usol] = ode23t(@(t,y) Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w), tspan, y0, opt);
    y0=usol(end,:);  
    i=endperiod/period;    
    NUM1=1024;    
    step=2*pi/omega/NUM1;
    tspan=[i*period: step: (i+1)*period];    
    Y1(mm)=norm(usol(:,2))/sqrt(1024);    
    [t,usol] = ode23t(@(t,y) Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w), tspan, y0, opt);
    Y1(mm)=norm(usol(:,2))/sqrt(1024);
    disp([' omega:', num2str(omega)]);
end

%%

%%
figure(2);
plot(f,Y1);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
%%

%%
function dy = Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w)

format long;

dy(1)=k1*(y(2)-y(1))-(knl_1*y(1)+knl_3*y(1)^3);
dy(2)=y(3);
dy(3)=(-c*y(3)-k1*(y(2)-y(1))+kAmp*cos(w*t))/m;
dy = dy';

end
%%

I want to do a frequency sweep for a DAE probelm. The code can work when I am not doing frequency sweep but calculating individual omega value. I know the reason why the code won't run is the new initial value didn't meet the requirement that the left side is 0 (it will be very close, but it may never reach).

dy(1)=k1*(y(2)-y(1))-(knl_1*y(1)+knl_3*y(1)^3);

Can anyone help me to solve the problem