def odeFun(t,y,**kwargs):

    """
    Contains system of differential equations.

    Arguments:
        t        : current time variable value
        y        : current state variable values (order matters)
        **kwargs : constant parameter values, interpolanting functions, etc.
    Returns:
        Dictionary containing dY/dt for the given state and parameter values
    """

    IW,IR = y # unpack state vars

    # Unpack variables passed through kwargs (see I thought I could avoid this
    # and I just made it messier)
    alp,bet,det,eps,N = kwargs['alp'],kwargs['bet'],kwargs['det'],kwargs['eps'], \
        kwargs['N']

    # ODEs
    dIW = ( bet * ( N + ( eps - 1 ) * IR - IW ) - det - alp ) * IW
    dIR = ( bet * ( N - ( eps + 1 ) * IW - IR ) - det ) * IR

    # Gather differential values in list (state variable order matters for
    # numerical solver)
    dy = [dIW,dIR]

    return dy