clear;c=1;k=10;m=1;fr=[0.5,1];
for n=1:2
    dydt=@(t,y) [y(2);(sin(fr(n)*2*pi*t)-c*y(2)-k*y(1))/m];
    [t,y]=ode45(dydt,[0,20],[0;0]);
    subplot(2,1,n);plot(t,y(:,1));grid('on');
    xlabel(['Frequency=',num2str(fr(n)),'[Hz]']);
    ylabel('Height[m]');ylim([-0.4,0.4]);
end
