あの程度の精度は出せるでしょ
出せないとしたら処理系の選択ミス

pi@raspi ~ $ cat simpson.py
from math import *

def f(x):
return x*(sin(x)+cos(x)+1.0)

N = 1000000

a = 0.0
b = 10.0**4.0
h = (b-a)/N

S = (h/3) * sum((f(h*i) + 4*f(h*(i+1)) + f(h*(i+2))) for i in range(0,N-1, 2))

print(S)
pi@raspi ~ $ python simpson.py
50006463.152