def simpson(x, f):
    """ 
    Compute a 1D definite integral using Simpson's rule.
    Parameters
    ----------
    f : function
        User defined function.
    x : numpy array
        Integration domain.
    Returns
    -------
    I : float
        Integration result.
    """   
    a = x[0]
    b = x[1]
    ya = f(a)
    yb = f((a+b)/2)
    yc = f(b)
    I = (b-a) * (ya + 4 * yb + yc) / 6
    return I