(* Algorithmus 7.2 *)
FUNCTION simpson(a,b,eps:real):real;
VAR h,s1,s2,s4,ta,tn : real; j,zh : integer;
    (* GLOBAL Funktion f: real *)
BEGIN
   h:=(b-a)/2; s1:=f(a)+f(b); s2:=0; s4:=f(a+h);
   tn:= h*(s1+4*s4)/3; zh:=2;
   REPEAT
     ta:=tn; zh:=2*zh; h:=h/2;
     s2:=s2+s4; s4:=0;
     j:=1;
     REPEAT
       s4:=s4+f(a+j*h); j:=j+2
     UNTIL j>zh;
     tn := h*(s1+2*s2+4*s4)/3;
     writeln(tf,zh:8,tn:25);
   UNTIL abs(ta-tn)<=eps*abs(tn);
   simpson := tn;
END;