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

@author: sam
"""

from scipy.special import p_roots
import numpy as np

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

def I(f,a,b,n):
    I = 0 #define I
    x,w = p_roots(n) #generate roots and weights
    for i in range(len(x)): #iterate through range of roots, 0 to n-1
        I += w[i]*f((b-a)/2*x[i]+(a+b)/2) #Gauss-Legendre quadrature formula
    I *= (b-a)/2 #multiplies I by (b-a)/2
    return I

print(I(f,-1,1,13)) #calls function I with function f, bounds [-1,1], and 13 roots since iteration goes up to 12