Oscillation with Damping and External Force

 Followings are the code that I wrote in Octave to creates all the plots shown in this page. You may copy these code and play with these codes. Change variables and try yourself until you get your own intuitive understanding.   < Code 1 >   clear;   t = 0:0.02:60; global m = 1.0; % 1.0 global k = 1.0; % 1.0 global c = 0.3; %0.3; global F0 = 0.50; %0.5; global w = 4.0; %2.5;    function dy_dt = f(y, t)   global m;   global k;   global c;   global F0;   global w;      dy_dt = zeros (2,1);   dy_dt(1) = y(2);   dy_dt(2) = -(c/m) .*y(2) -(k/m) .* y(1) + (F0/m) .* cos(w*t); endfunction   y2_init = 0.0; y1_init = 1.0;   y0 = [y1_init y2_init]; y = lsode ("f", y0, t);   hFig = figure(1,'Position',[300 300 780 300]);   subplot(1,5,[1 2 3]); plot(t,y(:,1),'r-','LineWidth',2,t,y(:,2),'b-','LineWidth',2); xlabel('time'); xlim([0 t(end)]); ylim([-2 2]); legend('y1(t)','y2(t)'); grid on;   subplot(1,5,[4 5]); plot(y(:,1),y(:,2),'r-','LineWidth',2);   xlabel('y(1)');ylabel('y(2)'); axis([-2 2 -2 2]); pbaspect([1 1 1]); tStr = sprintf("y1(0)=%0.2f,y2(0)=%0.2f \nm=%0.2f,k=%0.2f,c=%0.2f,F0=%0.2f,w=%0.2f",y1_init,y2_init,m,k,c,F0,w); title(tStr);