code: def get_coupled_ecc_eqns(y, t, mc, q):
      
      
      F = y[0]
      e = y[1]
      gamma = y[2]
      phase = y[3]
      
      #total mass
      m = (((1+q)**2)/q)**(3/5) * mc    
      
      dFdt = get_Fdot(F, mc, e)
      dedt = get_edot(F, mc, e)
      dgdt = get_gammadot(F, mc, q, e)
      dphasedt = 2*np.pi*F
       
      return np.array([dFdt, dedt, dgdt, dphasedt])
docs: Computes the coupled system of differential
      equations from Peters (1964) and Barack &
      Cutler (2004). This is a system of three variables:
      
      F: Orbital frequency [Hz]
      e: Orbital eccentricity
      gamma: Angle of precession of periastron [rad]
      phase0: Orbital phase [rad]
      
      :param y: Vector of input parameters [F, e, gamma]
      :param t: Time [s]
      :param mc: Chirp mass of binary [Solar Mass]
      :param q: Mass ratio of binary
      
      :returns: array of derivatives [dF/dt, de/dt, dgamma/dt, dphase/dt]