> restart:
> romberg:=proc(a,b,m,f)
> # [a,b] Integrationsintervall
> # m: Anzahl der Iterationen
> # f: Feld von Funktionswerten f(a),f(a+(b-a)/2^m),...,f(b)
> # option trace;  
> local  h,s,i,k,n,q,l,t;
> t:=array(1..m+1);
> h:=evalf(b-a); n:=1; l:=2^m;
> t[1]:=evalf(h*(f[1]+f[l+1])/2);
> for k from 1 to m do
>  s:=0; h:=h/2; n:=2*n; q:=1; l:=l/2;
>  for i from 1 by 2 to n-1 do s:=evalf(s+f[1+i*l]) od;
>  t[k+1]:=evalf(t[k]/2+s*h);
>  for i from k-1 by -1 to 0 do 
>    q:=4*q; t[i+1]:=evalf(t[i+2]+(t[i+2]-t[i+1])/(q-1))
>  od;
> od;
> evalf(t[1])
> end:
> a:=0: b:=1: m:=2: f:=array(1..5):
> for j from 0 to 4 do f[j+1]:=(a+j/4)^5 od:
> print(f);

                      [                 243    ]
                      [0, 1/1024, 1/32, ----, 1]
                      [                 1024   ]

> In:=romberg(a,b,m,f);

                          In := .1666666667

> Digits:=15:
> a:=0: b:=Pi/2: h:=b-a: m:=4: f:=array(1..17):
> for j from 0 to 16 do f[j+1]:=evalf(5*exp(2*(a+j*h/16))*cos(a+j*h/16)/(exp(Pi)-2)) od:
> print(f);

[.236510699382071, .286436274593345, .343535621355141,

    .407902361368191, .479247909030260, .556734605589705,

    .638760874395340, .722687951031825, .804496152493835,

    .878357002909815, .936105899063290, .966598479103435,

    .954932596766143, .881516999744180, .720967731790000,

    .440814276880832, 0]

> In:=romberg(a,b,m,f);

                        In := 1.00000000764186

> Digits:=15:
> a:=0: b:=Pi/2: h:=b-a: m:=6:
> f:=array(1..65):
> for j from 0 to 64 do f[j+1]:=evalf(5*exp(2*(a+j*h/64))*cos(a+j*h/64)/(exp(Pi)-2)) od:
> # print(f);
> In:=romberg(a,b,m,f);

                        In := 1.00000000000002

>