# -*- coding: utf-8 -*-
"""
Created on Sun Sep  4 14:18:14 2022

@author: Chris

Lab 6, Q7
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import p_roots

def f(x):
    f = np.cos(x)*np.exp(-2*x) 
    return f
    

def quadz(f,a,b,tol):
    
    NMAX = 100  #maximum number of points for quadrature
    Iold = 0
    
    for n in range(1,NMAX):
        # set up zeros and weights
        [xj,wj]=p_roots(n)
    
        I = (b-a)/2*sum(wj*f((b-a)/2*xj+(b+a)/2))
        #print(n,I)
        
        if n>3 and abs(I-Iold)<tol:
            return I
        
        Iold = I
        
    print('quadz: Quadrature did not converge to required tolerance.')
    return float("NaN")

result = quadz(f,1,3,1e-6)
print(result)