(*Compute the densities and nodes 
for Gaussian quadrature 
of degree (precision) 2  n + 
        1 over [a, b]; 
  generates n + 
        1 nodes and densities . 
Based on R . 
            Curto and L . 
            Fialkow
[
              Houston J .  Math .  17 (1991),  
              603 - 635], 
      MR 93  a : 
          47016*)         
Gaussian[
        n_, a_, b_] := (
For[i = 0, 
          i ≤ 
            2  n + 1, i++, 
bb[
              i]  =  (b^(i + 1) - a^(i + 1))/(i + 
                  1)]; 
For[i = 0, i ≤ n, 
          i = i + 1, 
For[j = 0, j ≤ n, 
            j = j + 1, 
Mtab[i, j] = 
              bb[i + j]]]; 
H = 
          Table[Mtab[r, s], {r, 0, n}, {s, 0, n}]; 
J  =  
          N[Inverse[H]]; 
For[i = 0,  i ≤ n, 
          i = i + 1, 
v[i]  =  
            bb[n + 1 + i]]; 
w  = 
          Table[v[i], {i, 0, n}]; 
z  =  
          J . w; 
p  =  0; 
For[i = 0, 
          i ≤ n, 
          i = i + 1, 
p  =  
            p  +  z[[i + 1]] t^i]; 
q  =  
          t^(n + 1) - p; 
roots  =  
          Solve[q == 0, t]; 
For[i = 1, 
          i ≤ n + 1, 
          i = i + 1, 
x[
              i - 1]  =  roots[[i]][[1]][[2]]]; 
For[ i = 0, i ≤ n, 
          i = i + 1, 
For[j = 0, j ≤ n, 
            j = j + 1, 
TTab[i, j] =  
              x[j]^i]]; 
V =  
          Table[TTab[rr, ss], {rr, 0, n}, {ss, 0, n}]; 
For[
          i = 0,  i ≤ n, 
          i = i + 1, 
vvv[i]  =  
            bb[i]]; 
w  = 
          Table[vvv[i], {i, 0, n}]; 
weights = 
          N[Inverse[V] . w];)

(*Gaussian quadrature estimate of
degree 2  
          n + 1 for integral of 
predefined function f[x]
over [a, b];  
      uses results of 
previous call to Gaussian*) 
G[
          j_] := (       (*j is ignored*) 
ss = 
            0; 
For[i = 0,  
            i ≤ n, i++, 
ss  =  
              ss  +  weights[[i + 1]] f[
                    x[i]]]; 
ss =  
            N[ss]);

a = 1

b = 3

n = 2

f[x_] := x^5

Gaussian[n, a, b]

G[1]

NIntegrate[f[x], {x, 1, 3}]

f[x_] :=  Sin[x]

G[1]

NIntegrate[f[x], {x, 1, 3}]

n = 4

Gaussian[n, a, b]

f[x_] :=  x^9

G[1]

NIntegrate[f[x], {x, 1, 3}]

f[x_] :=  Sin[5  x]

G[1]

NIntegrate[f[x], {x, 1, 3}]

Simpson[n_, a_, 
        b_] := (
NN = 
          2  n; 
h  =  (b - a)/
            NN; 
oddsum = 0; 
evensum = 
          0; 
For[i = 0,  i ≤ n - 1, 
          i = i + 1, 
oddsum  =  
            oddsum  +  
              f[a  +  (2  i + 1) 
                    h]; 
evensum  =   evensum  +  
              f[a + 2  i h]]; 
S =  
          h/3 ( 4 oddsum  +  2  evensum  -  f[a] + 
                f[b]))

N[Simpson[6, 1, 3]]

N[Simpson[7, 1, 3]]

n = 5

Gaussian[n, a, b]

G[1]

n = 6

Gaussian[n, a, b]

G[1]

n = 5

a = 0

b = 2

f[x_] := Exp[2  x] Sin[3  x]

Gaussian[n, a, b]

G[1]

NIntegrate[f[x], {x, 0, 2}]

N[Simpson[22, a, b]]

N[Simpson[23, a, b]]

Integrate[f[x], x]

FF[x_] :=  -(3/13) (E^(2 x)) Cos[3 x] + 
        (2/13) (E^(2 x)) Sin[3 x]

FF[x]

N[FF[2] - FF[0], 10]

Plot[f[x], {x, a, b}]