(*Implementacja metody Romberga.*)
(*fun - funkcja do scałkowania*)
(*eps - błąd, poniżej którego ma zejść program*)
(*lim1 i lim2 - granice całkowania*)

Romberg[fun_,eps_,{arg_,lim1_,lim2_}]:=Module[{h,calka,err,wynik},
h=lim2-lim1;
(*"zerowy rząd"*)
calka[0,0]=h/2*((fun/.arg->lim1)+(fun/.arg->lim2));
wynik=calka[0,0];
i=1;
err=calka[0,0];
(*Dałem arbitralnie 20, żeby program nie próbował pracować w nieskończoność.*)
While[i<20&&err>eps,
calka[i,0]=calka[i-1,0]/2+h*Sum[fun/.arg->(lim1+(2*k-1)*h/2^i),{k,1,2^(i-1)}]/2^i;
Do[calka[i,j]=calka[i,j-1]+(calka[i,j-1]-calka[i-1,j-1])/(4^j-1)
,{j,1,i}];
err=Abs[calka[i,i]-wynik];
wynik=calka[i,i];
i++;
];
(*Nie wyświetlam więcej cyfr po przecinku, niż wynika z zadanej wielkości eps*)
{N[wynik,Ceiling[-Log10[eps]]+RealDigits[N[wynik]][[2]]],i}
]

Romberg[Exp[x],10^-10,{x,0,1}]

N[E-1,10]