global omega_0 omega beta F0 m

m=1; %kg
beta=0.1; %wspolczynnik tlumienia
omega_0=10;% czestosc rezonansowa (rad/s)
F0=1; %amplituda sily wymuszajacej (N)
omegai=[9:0.05:11]; %zakres czestosci sily wymuszajace (rad/s)
t_zakres = [0, 150];  % przedzia³ czasowy
z0=[0,-3]; %warunki pocz±tkowe [x dla t=0,v dla t=0]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(omegai)
    omega=omegai(i);
[t, z] = ode23t('equation', t_zakres,z0); %wykonaj ode23
x = z(:,1); 
v= z(:,2); 
figure(1), plot(t,x), xlabel('t [s]', 'FontSize', 24, 'fontweight','bold'), ylabel('x [m]', 'FontSize', 24, 'fontweight','bold'),title(['\omega = ', num2str(omega),' rad/s']),set(gca, 'FontSize', 16)
pause(0.5)
A(i)=max(max(x(end-100:end)));
end
figure(2), plot(omegai,A), xlabel('\omega'), ylabel('A'), hold on
xlabel('\omega [rad/s]', 'FontSize', 24, 'fontweight','bold'), ylabel('A [m]', 'FontSize', 24, 'fontweight','bold'),set(gca, 'FontSize', 16), hold on


