code: def make_force_ES(q_vars, mass):
      
      # Build the potential energy fluxion; just one pair of bodies
      U =  U_ij(q_vars, mass, 0, 1)
      # Varname arrays for both the coordinate system and U
      vn_q = np.array([q.var_name for q in q_vars])
      vn_fl = np.array(sorted(U.var_names))
      # Permutation array for putting variables in q in the order expected by U (alphabetical)
      q2fl = np.array([np.argmax((vn_q == v)) for v in vn_fl])
      # Permutation array for putting results of U.diff() in order of q_vars
      fl2q = np.array([np.argmax((vn_fl == v)) for v in vn_q])
      # Return a force function from this potential
      force_func = lambda q: -U.diff(q[q2fl]).squeeze()[fl2q]
      return force_func
docs: Fluxion with the potential energy of the earth-sun sytem