PROGRAM ROMBERG(INPUT, OUTPUT);                   
(* NUMERISCHE INTEGRATION DER FUNKTION F.       *)
(* ES WIRD EINE SUMMIERTE SEHNEN-TRAPEZFORMEL   *)
(* MIT 2^J STUETZSTELLEN BERECHNET. ZUSAETZLICH *)
(* WIRD EINE LINEARKOMBINATION DER ZUGEHOERIGEN *)
(* NAEHERUNGEN BENUTZT, UM FEHLERTERME HOEHERER *)
(* ORDNUNG ZU ELIMINIEREN.                      *)
(* QUELLE: STUMMEL/HAINER: PRAKTISCHE MATHEMATIK*)
(*         TEUBNER STUDIENBUECHER MATHEMATIK    *)
                                                  
(*$R+ INDEXGRENZEN PRUEFEN                      *)
 CONST MX=7;             (* J<=MX               *)
       LEN=13;                                    
 VAR A,B  : REAL;        (* INTEGRATIONSGRENZEN *)
     S0,F4,W: REAL;                               
     H    : REAL;        (* ABSTAND DER STUETZS.*)
     NJ   : INTEGER;     (* NJ=2^J STUETZSTELLEN*)
     I,J,K: INTEGER;     (* LAUFVARIABLEN       *)
     S    : ARRAY[0..MX] OF REAL;                 
     CH   : CHAR;                                 
                                                  
FUNCTION F(X:REAL):REAL;                          
BEGIN                                             
 F:=X*X*EXP(X)                                    
END; (* F *)                                      
                                                  
BEGIN                                             
  WRITELN("INTEGRATIONSGRENZEN:");                
  WRITE("VON:");READLN(A);                        
  WRITE("BIS:");READLN(B);                        
  WRITE("AUSGABE AUF DEN DRUCKER? ");READLN(CH);  
  IF CH IN ["J","Y"] THEN OPEN(OUTPUT,4,0);       
  NJ:=1;                                          
  FOR J:= 0 TO MX DO                              
    BEGIN                                         
      S0:=0.0;                                    
      H:= (B-A) / NJ;                             
      FOR I:= 0 TO NJ DO                          
        BEGIN                                     
          W:= A+I*H;            (* STUETZSTELLE *)
          IF (I=0) OR (I=NJ) THEN                 
             S0:= S0+F(W)/2                       
          ELSE                                    
             S0:= S0+F(W)                         
        END;                                      
      S[J]:= S0*H; F4:=4;       (* F4 = 4^K+1   *)
      FOR K:= J-1 DOWNTO 0 DO                     
        BEGIN                                     
          S[K]:=(F4*S[K+1]-S[K])/(F4-1);          
          F4:=F4*4                                
        END;                                      
      (* DRUCKE ZWISCHENSUMMEN *)                 
      WRITE(J:3);                                 
      FOR K:=J DOWNTO 0 DO                        
        WRITE("Ý",S[K]:LEN);                      
      WRITELN;                                    
      NJ:=NJ*2                                    
    END; (* FOR J ...*)                           
  CLOSE(OUTPUT)                                   
END.