code: def leapfrog(func,yo,t,args=(),rtol=1.49012e-12,atol=1.49012e-12):
      
      #Initialize
      qo= yo[0:len(yo)//2]
      po= yo[len(yo)//2:len(yo)]
      out= nu.zeros((len(t),len(yo)))
      out[0,:]= yo
      #Estimate necessary step size
      dt= t[1]-t[0] #assumes that the steps are equally spaced
      init_dt= dt
      dt= _leapfrog_estimate_step(func,qo,po,dt,t[0],args,rtol,atol)
      ndt= int(init_dt/dt)
      #Integrate
      to= t[0]
      for ii in range(1,len(t)):
          for jj in range(ndt): #loop over number of sub-intervals
              #This could be made faster by combining the drifts
              #drift
              q12= leapfrog_leapq(qo,po,dt/2.)
              #kick
              force= func(q12,*args,t=to+dt/2)
              po= leapfrog_leapp(po,dt,force)
              #drift
              qo= leapfrog_leapq(q12,po,dt/2.)
              #Get ready for next
              to+= dt
          out[ii,0:len(yo)//2]= qo
          out[ii,len(yo)//2:len(yo)]= po
      return out
docs: NAME:
         leapfrog
      PURPOSE:
         leapfrog integrate an ode
      INPUT:
         func - force function of (y,*args)
         yo - initial condition [q,p]
         t - set of times at which one wants the result
         rtol, atol
      OUTPUT:
         y : array, shape (len(y0), len(t))
         Array containing the value of y for each desired time in t, \
         with the initial value y0 in the first row.
      HISTORY:
         2011-02-02 - Written - Bovy (NYU)