function simpson(f, a, b, level, level_max, precision)
    level += 1
    h = b - a
    c = (a + b) / 2
    one_s = h * (f(a) + 4. * f(c) + f(b)) / 6.
    d = (a + c) / 2
    e = (c + b) / 2
    two_s = h * (f(a) + 4. * f(d) + 2. * f(c) + 4. * f(e) + f(b)) / 12.
    if level >= level_max
        simpson_result = two_s
        print("Max level reached")
        return simpson_result
    else
        if abs(two_s - one_s) < 15. * precision
            return two_s + (two_s - one_s) / 15.
        else
            left_s = simpson(f, a, c, level, level_max, precision / 2.)
            right_s = simpson(f, c, b, level, level_max, precision / 2.)
            print(left_s)
            return left_s + right_s
        end
    end
end

function simptest(x)
    fn = x^2
    return fn
end

simpson(simptest, -2, 2, 0.000000000000001, 1, 10) #returns 5.333333333333.