function no_result = run_resonance()
   
   global omega A g L
   
   %initial conditions
   x0=0.0;
   v0=0.0;
   t0=0.0;
   
   %timing
   dt = 0.2
   duration = 100
   N_steps = duration/dt
   
   % forcing function characteristics
   omega = 1.0;     % forcing function frequency (sec^{-1})
   A = 0.1;        %  forcing function amplitude (sec^{-2})
   
   % Note: A=0.01 gives a typical Lorentzian resonance profile, typical
   %        of the linear damped/driven oscillator. A=0.1 gives a very
   %        asymmetric profile due to the nonlinearity of the pendulum.
   
   N_omegas = 50;
   domega = 0.01;
   omega_0 = 0.7;
   omega_array = omega_0 + (0:N_omegas-1)*domega;
   for i = 1:N_omegas
      omega = omega_array(i);
      y = rk2([x0,v0], t0, dt, N_steps, 'pendulum2');
      theta = y(:,1);
      theta_dot = y(:,2);
      energy = g*L*(1 - cos(theta)) + 0.5*(L*theta_dot).^2;
      energy_array(i) = max(energy);
   end
   
   time = t0:dt:duration;
   hold off
   plot(omega_array, energy_array,'r')
