function int=simpson_quad(fa, fab2, fb, a, b)   
    int = (b - a) / 6.0 * (fa + 4*fab2 + fb);   
end