def custom_splitting_operator(fa, fb, ai, bi):
    """
    Produces a splitting rule (evolution operator) from the given evolution operators and coefficients
    fa: first evolution operator
    fb: second evolution operator
    ai: collection of float
    bi: collection of float

    Preconditions:
        -length of ai and bi match (add zeros yourself, if necessary)
        -sum of values in ai/bi are approximately 1
    """
    
    if abs(sum(ai)-1)>0.0001 or abs(sum(bi)-1)>0.0001:
        raise Exception("Coefficients do not add up to 1...")
    
    def operator(y0, dt):
        result = y0
        count = 0
        for a, b in zip(ai, bi):
            result = fb(fa(result, a*dt), b*dt)
        return result

    return operator