> restart;
> # ADAPTIVE QUADRATURE ALGORITM  4.3 #    
> # To approximate I = integral ( ( f(x) dx ) ) from a to b to within
> # a given tolerance TOL:
> #
> # INPUT:    endpoints a, b; tolerance TOL; limit N to number of levels
> #
> # OUTPUT:   approximation APP or message that N is exceeded.
> alg043 := proc() local F, OK, AA, BB, EPS, N, CNT, APP, I, TOL, A, H, FA, FC, FB, S, L, FD, FE, S1, S2, V, LEV;
> printf(`This is Adaptive Quadrature with Simpson's Method.\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`);
> AA := scanf(`%f`)[1];
> BB := scanf(`%f`)[1];
> if AA > BB then
> printf(`Lower limit must be less than upper limit\n`);
> else
> OK := TRUE;
> fi;
> od; 
> OK := FALSE;
> while OK = FALSE do
> printf(`Input tolerance.\n`);
> EPS := scanf(`%f`)[1];
> if EPS > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be positive.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of levels.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> OK := TRUE;
> else
> printf(`Number must be positive\n`);
> fi;
> od;
> if OK = TRUE then
> CNT := 0;
> OK := TRUE;
> # Step 1
> APP := 0;
> I := 1;
> TOL[I] := 10*EPS;
> A[I] := AA;
> H[I] := 0.5*(BB-AA);
> FA[I] := F(AA);
> CNT := CNT+1;
> FC[I] := F((AA+H[I]));
> CNT := CNT+1;
> FB[I] := F(BB);
> CNT := CNT+1;
> # Approximation from Simpson's method for entire interval
> S[I] := H[I]*(FA[I]+4*FC[I]+FB[I])/3;
> L[I] := 1;
> # Step 2
> while I > 0 and OK = TRUE do
> # Step 3
> FD := F((A[I]+0.5*H[I]));
> CNT := CNT+1;
> FE := F((A[I]+1.5*H[I]));
> CNT := CNT+1;
> # Approximations from Simpson's method for halves of intervals
> S1 := H[I]*(FA[I]+4*FD+FC[I])/6;
> S2 := H[I]*(FC[I]+4*FE+FB[I])/6;
> # Save data at this level
> V[0] := A[I];
> V[1] := FA[I];
> V[2] := FC[I];
> V[3] := FB[I];
> V[4] := H[I];
> V[5] := TOL[I];
> V[6] := S[I];
> LEV := L[I];
> # Step 4
> # Delete the level
> I := I-1;
> # Step 5
> if abs(S1+S2-V[6]) < V[5] then
> APP := APP+(S1+S2);
> else
> if LEV >= N then
> OK := FALSE;
> # Procedure fails
> else
> # Add one level
> # Data for right half subinterval
> I := I+1;
> A[I] := V[0]+V[4];
> FA[I] := V[2];
> FC[I] := FE;
> FB[I] := V[3];
> H[I] := 0.5*V[4];
> TOL[I] := 0.5*V[5];
> S[I] := S2;
> L[I] := LEV+1;
> I := I+1;
> # Data for left half subinterval
> A[I] := V[0];
> FA[I] := V[1];
> FC[I] := FD;
> FB[I] := V[2];
> H[I] := H[I-1];
> TOL[I] := TOL[I-1];
> S[I] := S1;
> L[I] := L[I-1];
> fi;
> fi;
> od;
> if OK = FALSE then
> printf(`Level exceeded. Method failed to produce accurate approximation.\n`);
> else
> printf(`\nThe integral of F from %12.8f to %12.8f is\n`, AA, BB);
> printf(`%12.8f to within %14.8e\n`, APP, EPS);
> printf(`The number of function evaluations is: %d\n`, CNT); 
> fi;
> fi;
> RETURN(0);
> end;
> alg043();