# -*- coding: utf-8 -*-
"""
Created on Sun Sep  4 17:35:57 2022

@author: Chris
Lab 6, Q3
"""

import numpy as np
import matplotlib.pyplot as plt

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


def qtrapn(f,a,b,N):
    # perform trapezoidal rule on a 1D function f, 
    # on the interval [a,b], with N subintervals
    
    h = (b-a)/N
    
    xj = np.array([a+h*j for j in range(N+1)])

    I = h/2*(f(xj[0]) + 2*sum(f(xj[1:-1])) + f(xj[-1]))
    
    return I

def qsimpn(f,a,b,N):
    # perform simpson's composite rule on a 1D function f, 
    # on the interval [a,b], with N subintervals
    
    h = (b-a)/N
    
    xj = np.array([a+h*j for j in range(N+1)])

    I = h/3*(f(xj[0]) + 4*sum(f(xj[1:-1:2]))
                      + 2*sum(f(xj[2:-1:2]))
                      + f(xj[-1]))
    
    return I

a = 0
b = 2

for p in range(2,11):
    n = 2**p

    result = qsimpn(f,a,b,n)

    print(n,result)