code: def direct_nbody(q,p,m,t,pot=None,softening_model='plummer',
                   softening_length=None,
                   atol=None,rtol=None):
      
      #Set up everything
      if softening_model.lower() == 'plummer':
          softening= _plummer_soft
      out= []
      out.append([q,p])
      qo= q
      po= p
      #Determine appropriate stepsize
      dt= t[1]-t[0]
      ndt= 1
      to= t[0]
      #determine appropriate softening length if not given
      if softening_length is None:
          softening_length= 0.01
      #Run simulation
      for ii in range(1,len(t)):
          print(ii)
          for jj in range(ndt): #loop over number of sub-intervals
              (qo,po)= _direct_nbody_step(qo,po,m,to,dt,pot,
                                          softening,(softening_length,))
              #print(qo)
              to+= dt
          out.append([qo,po])
      #Return output
      return out
docs: NAME:
         direct_nbody
      PURPOSE:
         N-body code using direct summation for force evaluation
      INPUT:
         q - list of initial positions (numpy.ndarrays)
         p - list of initial momenta (numpy.ndarrays)
         m - list of masses
         t - times at which output is desired
         pot= external potential (galpy.potential or list of galpy.potentials)
         softening_model=  type of softening to use ('plummer')
         softening_length= (optional)
      OUTPUT:
         list of [q,p] at times t
      HISTORY:
         2011-02-03 - Written - Bovy (NYU)