> restart;
> # ROMBERG ALGORITHM 4.2
> #     
> # To approximate I = integral ( f(x) dx ) from a to b:
> #
> # INPUT:   endpoints a, b; integer n.
> #
> # OUTPUT:  an array R. ( R(2,n) is the approximation to I. )
> #
> # R is computed by rows; only 2 rows saved in storage   
> alg042 := proc() local F, OK, A, B, N, H, R, I, SUM, M, K, J, L;
> printf(`This is Romberg Integration.\n\n`);
> printf(`Input the function F(x) in terms of x\n`);
> printf(`For example: cos(x)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input lower limit of integration and `);
> printf(`upper limit of integration\n`);
> printf(`separated by a blank\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A > B then
> printf(`Lower limit must be less than upper limit\n`);
> else
> OK := TRUE;
> fi;
> od; 
> OK := FALSE;
> while OK = FALSE do
> printf(`Input number of rows - no decimal point.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> OK := TRUE;
> else
> printf(`Number must be a positive integer\n`);
> fi;
> od;
> if OK = TRUE then
> # STEP 1
> H := B-A;
> R[0,0] := (F(A)+F(B))/2*H;
> # STEP 2
> printf(`Initial Data:\n`);
> printf(`Limits of integration = [%12.8f, %12.8f]\n`, A, B);
> printf(`Number of rows = %3d\n`, N);
> printf(`\nRomberg Integration Table:\n`);
> printf(`\n%12.8f\n\n`, R[0,0]);
> # STEP 3
> for I from 2 to N do
> # STEP 4
> # approximation from Trapezoidal method
> SUM := 0;
> M := 2^(I-2);
> for K from 1 to M do
> SUM := SUM+F(A+(K-0.5)*H);
> od;
> R[1,0] := (R[0,0]+H*SUM)/2;
> # STEP 5
> # extrapolation
> for J from 2 to I do
> L := 2^(2*(J-1));
> R[1,J-1] := R[1,J-2]+(R[1,J-2]-R[0,J-2])/(L-1);
> od;
> # STEP 6
> for K from 1 to I do
> printf(` %11.8f`,R[1,K-1]);
> od;
> printf(`\n\n`);
> # STEP 7
> H := H/2;
> # since only two rows are kept in storage, this step is 
> # to prepare for the next row
> # update row 1 of R
> # STEP 8
> for J from 1 to I do
> R[0,J-1] := R[1,J-1];
> od;
> od;
> fi;
> RETURN(0);
> end;
> alg042();