<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def find_discrete(start_time, end_time, f, epsilon=EPSILON, num=12):
 """Find the times when a function changes value. Searches between ``start_time`` and ``end_time``, which should both be :class:`~skyfield.timelib.Time` objects, for the occasions where the function ``f`` changes from one value to another. Use this to search for events like sunrise or moon phases. A tuple of two arrays is returned. The first array gives the times at which the input function changes, and the second array specifies the new value of the function at each corresponding time. This is an expensive operation as it needs to repeatedly call the function to narrow down the times that it changes. It continues searching until it knows each time to at least an accuracy of ``epsilon`` Julian days. At each step, it creates an array of ``num`` new points between the lower and upper bound that it has established for each transition. These two values can be changed to tune the behavior of the search. """

    ts = start_time.ts
    jd0 = start_time.tt
    jd1 = end_time.tt
    if jd0 >= jd1:
        raise ValueError('your start_time {0} is later than your end_time {1}'
                         .format(start_time, end_time))

    periods = (jd1 - jd0) / f.rough_period
    if periods < 1.0:
        periods = 1.0

    jd = linspace(jd0, jd1, periods * num // 1.0)

    end_mask = linspace(0.0, 1.0, num)
    start_mask = end_mask[::-1]
    o = multiply.outer

    while True:
        t = ts.tt_jd(jd)
        y = f(t)

        indices = flatnonzero(diff(y))
        if not len(indices):
            return indices, y[0:0]

        starts = jd.take(indices)
        ends = jd.take(indices + 1)

        # Since we start with equal intervals, they all should fall
        # below epsilon at around the same time; so for efficiency we
        # only test the first pair.
        if ends[0] - starts[0] <= epsilon:
            break

        jd = o(starts, start_mask).flatten() + o(ends, end_mask).flatten()

    return ts.tt_jd(ends), y.take(indices + 1)