#!/usr/bin/python3
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import p_roots #for pulling zeroes of Legendre orthogonal polynomials


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

def quadz(f,a,b,tol):
    maxn=100
    I0=0
    for n in range(1,maxn):
        x,w=p_roots(n)
        I=(b-a)/2*sum(w*f((b-a)/2*x+(a+b)/2))
        print(n,I)
        if abs(I0-I)<tol:
            return I
        I0=I
    print("tol unreached")
    return None






quadz(f,0,2,1e-8)