Trying to program a chebyshev quadrature solution of the integral of e^(-x^2)between 0 and 1 using python

I can't figure out what I'm doing wrong. I'm ending up with an answer of 0.84965 while it should be much closer to 0.746824

import numpy as np
import math

def f(x):
    return np.e**(-x**2)

def chebyshev(a,b,n):
    h=(b-a)/n
    s=0
    r=(b-a)/2
    for i in range(n+1):
        l=i*h
        if i>n/2:
            l=h*(n/2-(i-n/2))
        val=np.cos(np.arcsin(l/r))*r
        s+=f(val)

    return s*h
print(chebyshev(0,1,1000))

If I use the equation from page 11 of these notes for the Chebyshev-Gauss approximation, and do:

from math import pi, cos, sqrt, exp


def chebyshev_node(i, n):
    return cos(0.5 * pi * (2 * i - 1) / n)


def chebyshev_solver(f, a, b, n):
    d = (b - a)
    c = 0.5 * pi * d / n

    s = 0.0
    for i in range(1, n + 1):
        cn = chebyshev_node(i, n)
        v = 0.5 * (cn + 1) * d + a
        s += f(v) * sqrt(1 - cn**2)

    return s * c


def efunc(x):
    return exp(-(x**2))


print(chebyshev_solver(efunc, 0, 1, 1000))

it gives 0.7468244140713791, which seems to match your expected solution.

Update: Just to note, the above could all be vectorised with NumPy as:

import numpy as np


def chebyshev_nodes(n):
    return np.cos(0.5 * np.pi * (2 * np.arange(1, n + 1) - 1) / n)


def chebyshev_solver(f, a, b, n):
    d = (b - a)
    c = 0.5 * np.pi * d / n

    cn = chebyshev_nodes(n)
    v = 0.5 * (cn + 1) * d + a
    return c * np.sum(f(v) * np.sqrt(1 - cn**2))


def efunc(x):
    return np.exp(-(x**2))

print(chebyshev_solver(efunc, 0, 1, 1000))