expr[x_]:=a*(x^2)+b*x+c;
eqns={fMinus==expr[x0-Δx],
f0==expr[x0],fPlus==expr[x0+Δx]};
solabc=Solve[eqns,{a,b,c}];
integrand=expr[x]/.solabc[[1]];
Integrate[integrand,{x,x0-Δx,x0+Δx}]

3-(((-1))^(#))&/@Range[11]

weights[n_]:=Append[Prepend[3-(((-1))^(#))&/@Range[n],1],1]
ClearAll[simpsonInt]
simpsonInt[f_,{x1_,x2_},n_?OddQ]:=Module[{w,nOfValues,Δx,xValues,functionValues},
w=weights[n];
nOfValues=n+2;
Δx=((x2-x1)/(nOfValues-1.));
xValues=Range[x1,x2,Δx];
functionValues=f/@xValues;
(functionValues.w)*((Δx)/3.)
]
simpsonInt[Sin,{0,0.5Pi},15]-1

ClearAll[simpsonInt]
simpsonInt[f_,{x1_,x2_},n_?OddQ]:=Module[{w,nOfValues,Δx,xValues,functionValues},
w=3&/@Range[n+1];
nOfValues=n+2;
Δx=((x2-x1)/(nOfValues-1.));
xValues=Drop[Range[x1,x2,Δx],-1];
functionValues=f/@xValues;
(functionValues.w)*((Δx)/3.)
]
simpsonInt[Sin,{0,0.5Pi},15111]-1

Plot[Evaluate@Table[PDF[NormalDistribution[0,σ],x],{σ,{.75,1,2}}],{x,-6,6},Filling->Axis]