The numerical integration follows the exaplanation on the online notes here:

https://engcourses-uofa.ca/books/numericalanalysis/numerical-integration/introduction/

The code used:

import sympy as sp
from scipy.integrate import quad
import matplotlib.pyplot as plt
import numpy as np
import base64
import io
import math


#This function is used to create an image from the plot created using matplotlib.plt
def plt_show(plt, width=500, dpi=100):
    # Converts matplotlib plt to data string
    # dpi (dots per inch) is the resolution of the image
    # width is width of image in pixels
    bytes = io.BytesIO()
    plt.savefig(bytes, format='png', dpi=dpi)  # Save as png image
    plt.close()
    bytes.seek(0)
    base64_string = "data:image/png;base64," + \
        base64.b64encode(bytes.getvalue()).decode("utf-8")
    return ""


#The following functions calculate the integrals and the coordinates for plotting the visual curves representing the integration
def I1(f, a, b, n):
  h = (b - a)/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x_coor.extend([(a+i*h),(a+(i+1)*h)])
    y_coor.extend([f(a+i*h),f(a+(i)*h)])
    i=i+1
  return [sum([f(a + i*h)*h for i in range(int(n))]),x_coor,y_coor]
def I2(f, a, b, n):
  h = (b - a)/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x_coor.extend([(a+i*h),(a+(i+1)*h)])
    y_coor.extend([f(a+(i+0.5)*h),f(a+(i+0.5)*h)])
    i=i+1
  return [sum([f(a + (i + 1/2)*h)*h for i in range(int(n))]),x_coor,y_coor]
def I3(f, a, b, n):
  h = (b - a)/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x_coor.extend([(a+i*h),(a+(i+1)*h)])
    y_coor.extend([f(a+(i+1)*h),f(a+(i+1)*h)])
    i=i+1
  return [sum([f(a + (i + 1)*h)*h for i in range(int(n))]),x_coor,y_coor]
def IT(f, a, b, n):
  h = (b - a)/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x_coor.extend([(a+i*h),(a+(i+1)*h)])
    y_coor.extend([f(a+(i)*h),f(a+(i+1)*h)])
    i=i+1
  return [sum([(f(a + i*h) + f(a + (i + 1)*h))*h/2 for i in range(int(n))]),x_coor,y_coor]
def IS1(f, a, b, n):
  h = (b - a)/2/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x1=np.linspace((a+2*i*h),(a+(i+1)*2*h))
    xi1=a+i*2*h
    xi=a+(i+1)*2*h
    xmi=a+(i+0.5)*2*h
    y1=f(xi1)*(x1-xmi)*(x1-xi)/2/h**2-f(xmi)*(x1-xi1)*(x1-xi)/h**2+f(xi)*(x1-xi1)*(x1-xmi)/2/h**2
    x_coor.extend(x1)
    y_coor.extend(y1)
    i=i+1
  return [h/3*sum([(f(a + (2*i)*h) + 4*f(a + (2*i + 1)*h) + f(a + (2*i + 2)*h)) for i in range(int(n))]),x_coor,y_coor]
def IS2(f, a, b, n):
  h = (b - a)/3/n
  x_coor=[]
  y_coor=[] 
  i=0;
  while i< int(n):
    x1=np.linspace((a+i*3*h),(a+(i+1)*3*h))
    xi1=a+i*3*h
    xli=xi1+h
    xri=xli+h
    xi=xri+h
    #y1=x1**2+2.0
    #print(y1)
    y1=f(xi1)*(x1-xli)*(x1-xri)*(x1-xi)/-6/h**3+f(xli)*(x1-xi1)*(x1-xri)*(x1-xi)/2/h**3-f(xri)*(x1-xi1)*(x1-xli)*(x1-xi)/2/h**3+f(xi)*(x1-xi1)*(x1-xli)*(x1-xri)/6/h**3
    x_coor.extend(x1)
    y_coor.extend(y1)
    i=i+1
  return [3*h/8*sum([(f(a + (3*i)*h) + 3*f(a + (3*i + 1)*h) + 3*f(a + (3*i + 2)*h) + f(a + (3*i + 3)*h)) for i in range(int(n))]),x_coor,y_coor]


#The following function is used to select the appropriate case
def switch(Method,f,a,b,n):
  tester={
"Rectangle Method 1 (Left Corner)":I1,
"Rectangle Method 2 (Midpoint)":I2,
"Rectangle Method 3 (Right Corner)":I3,
"Trapezoid Rule":IT,
"Simpson's One Third Rule":IS1,
"Simpson's Three Eightth Rule":IS2,
  }
  return tester.get(Method,"none found")(f,a,b,n)


#The following function is used for rounding to significant digits. 
def rounds(a_number,significant_digits):
  small_value=1.0e-9
  try:
    rounded_number =  round(a_number, significant_digits - int(math.floor(math.log10(max(abs(a_number),abs(small_value))))) - 1)
  except:
    rounded_number = "N/A"
  return rounded_number



#This is the main function
def main(inputs):
    f=inputs['func']
    a=inputs['a']
    b=inputs['b']
    Method=inputs['method']
    n=inputs['n']
    Err=False
    img=False
    Err2=False
    try:
        g=sp.sympify(f)
    except:
        Err=True
    #Here I use sympy plot which is has matplotlib as its backend.
    if not Err:
        x=sp.symbols('x')
        try:
            sp.plot(g,(x,a,b))
            img=plt_show(plt,500)
            f=sp.lambdify(x,g)
        except:
            Err2=True
    if not Err and not Err2:
        res,err=quad(f,a,b)
        #Exact="The numerical result is {:f} (+-{:g})".format(res, err))
        #print(switch(Method,f,a,b,n))
        Integration,x_coor,y_coor=switch(Method,f,a,b,n)
        percent=rounds((Integration-res)/res*100,3)
        plt.fill_between(x_coor,0,y_coor,color='black',facecolor='orange')
        x1=np.linspace(a,b)
        plt.plot(x1,f(x1))
        curve=plt_show(plt,500)
    else:
        print("error")
    return {"plot":img, "equation":g,"Error":Err,"ErrorPlot":Err2,"a":a,"b":b,"n":n,"Method":Method,"Integration":rounds(Integration,5),"Res":res,"Err":err,"Percent":percent,"Curve":curve}

{% if outputs.Error %}
Please Enter a valid expression in x
{% elif outputs.ErrorPlot %}
Function could not be plotted. Please enter a valid expression in variable x
{% else %}
The following figure provides the plot for the function f(x)={{outputs.equation}} between {{outputs.a}} and {{outputs.b}}
{{outputs.plot|safe}}
The python scipy package provides an accurate estimate of the ingral which is {{outputs.Res}}+/-{{outputs.Err}}
Using {{outputs.Method}}, numerical integration with {{outputs.n}} intervals yields {{outputs.Integration}}, which is {{outputs.Percent}}% different from the scipy result.
The following curve provides the plot of the function and the area calculated according to {{outputs.Method}}:
{{outputs.Curve|safe}}
{% endif %}

The output page:

{% if outputs.Error %}
Please Enter a valid expression in x
{% elif outputs.ErrorPlot %}
Function could not be plotted. Please enter a valid expression in variable x
{% else %}
The following figure provides the plot for the function f(x)={{outputs.equation}} between {{outputs.a}} and {{outputs.b}}
{{outputs.plot|safe}}
The python scipy package provides an accurate estimate of the ingral which is {{outputs.Res}}+/-{{outputs.Err}}
Using {{outputs.Method}}, numerical integration with {{outputs.n}} intervals yields {{outputs.Integration}}, which is {{outputs.Percent}}% different from the scipy result.
The following curve provides the plot of the function and the area calculated according to {{outputs.Method}}:
{{outputs.Curve|safe}}
{% endif %}