#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 10 14:37:07 2024

@author: sam
"""

def f(x):
    return x**3 + 1

def qtrapn(f,a,b,N):
    h = (b-a)/N
    I = f(a) + f(b)
    for i in range(1,int(N)):
        xi = a + i*h
        I += 2*f(xi)
    I *= h/2
    return I

def qtrapz(f,a,b,tol): #this doesnt work
    N = 2
    In = 0
    correction = 0
    h = (b-a)/N
    I = qtrapn(f,a,b,N)
    while abs(In - 6) >= tol:
        for i in range(1,int(N)):
            xi = a + i*h
            xij = (xi + i*h)/2
            correction += f(xij)
        In = I/2 + (h/2)*correction
        xi += xij
        N *= 2
        h *= 1/2
    return In

result = qtrapz(f,0,2,1e-3)
print(result)