function dydt = oscillator_eq (t,y)
     g = 1;
     R = 1;

     % This shows how to do a second order ODE
     % This represents a system of equations
     %   y(1)' = y(2)
     %   y(2)' = (2*cos(y(1))-(g/R))*sin(y(1))
     % Here y(1) = theta, so the first equation says y(2) = theta'
     % and the second equation says 
     %   y(2)' = theta" = (2*cos(theta)-(g/R))*sin(theta)

     %     y(1)'         y(2)'
   dydt = [y(2); (2*cos(y(1))-(g/R))*sin(y(1))];
   
   % Here is a third order example:
   % The equation y''' + 2y' = 3y would be coded:
   %    dydt = [y(2); y(3); 3*y(1) - 2*y(2)];
   % because, using y = y(1)
   %    y(3)' = 3*y(1) - 2*y(2)
   %    y(2)' = y(3)
   %    y(1)' = y(2)
    
